A»«730id 


mm 


CS  104 


A  NUMERICAL  INVESTIGATION 
OF  THE 

SIMPLEX  METHOD 


RICHARD  H.  BARTELS 


TECHNICAL  REPORT  NO.  CS  104 
JULY  31,  1968 


AUG  1  a  U68 


COMPUTER  SCIENCE  DEPARTMENT  r 
School  of  Humanities  and  Sciences 
STANFORD  UNIVERSITY 


Reproduced  by  Ihe 

CLEARINGHOUSE 
for  Federal  Scientific  &  Technical 
Informefion  Springfield  Ve.  22151 


'*/ 


1 


A  Numerical  Investigation 
of  the 

Simplex  Method* 


By 


Richard  H.  Bartels 


Computer  Science  Department 


Stanford  University 


^Research  sponsored  by  the  National  Science  Foundation  under  contract 
NSF  GP  5962  and  by  the  Office  of  Naval  Research  under  contract 
NOOOlU -67 -A-0112-0029 . 


Acknowledgments 

I  am  happy  at  having  this  chance  to  express  my  gratitude  to: 

who,  first  as  my  fiancee  and  then  as  my  wife,  lifted  my 
spirits  through  the  final  agonies  of  doctoral  candidacy; 

Qene  Golub,  my  thesis  advisor  and  my  friend,  who  has  given  me 
measureless  help,  academically  and  privately,  throughout  my  stay  at 
Stanford ; 

Professor  George  E.  Forsythe,  whom  I  thank  for  his  aid  and  advice 
to  me  during  the  course  of  my  academic  progress; 

Professor  George  B.  Dantzig,  whose  words  of  encouragement  to  me 
were  very  welcome; 

And  Miss  Suzanne  Stone,  who  under  the  pall  of  a  deadline,  did  all 
of  the  typing  which  I  required  and  maintained  her  composure  throughout. 

I  would  like  to  note  that  my  work  at  the  Computer  Science  Depart¬ 
ment  has  been  financially  supported  by  the  National  Science  Foundation 
and  the  Office  of  Naval  Research. 

And  I  give  my  thanks  to  Doctor  James  H.  Wilkinson  and,  again,  to 
George  B.  Dantzig.  They  laid  the  foundations. 


iii 


Table  of  Content* 


Chapter  Page 

1.  Introduction  and  Summary  of  Results  .  1 

2.  Review  of  the  Simplex  Method  of  Linear  Programming  •  •  .  8 

3.  Standard  Computer  Implementation*  of  the  Simplex  Method  .  13 

4.  Round-off  Error  Analyaia  of  Gauaa-Jordan  Elimination  .  •  18 

3*  Error  Bounda  for  Linear-ayatem  Solutions  Found  by 

Explicit  Inverse  .  23 

6.  The  Modified  LU  Implementation  of  the  Simplex  Method  .  .  34 

7.  Round-off  Error  Analysis  of  the  Modified  LU 

Implementation . • . 4o 

8.  The  Standard  LU  Implementation  of  the  Simplex  Method  •  •  75 

9.  Error  Bounds  for  Linear-system  Solutions  Found  by 

LU  Decomposition . 84 

10.  Employment  of  Error  Bounds  in  the  Simplex-method  Cycle  •  93 

11.  Comparison  of  Implementations  by  Operation  Counts  ....  103 

Bibliography  .  118 


V 


I  »  I  * 


-  r  wiraiiww  •?* 


1.  Introduction  and  Summary  of  Beaulta 

Since  19^7  linear  progranalng  problema  have  been  growing  ever  more 
interesting  to  business,  industry,  government,  and  the  military.  Today 
few  facets  of  the  economy  remain  uninfluenced  by  decisions  made  with 
the  aid  of  linear  programming.  The  increasing  practical  use  of  linear 
programming  since  19^7  is  linked  with  the  rise  of  electronic  computers, 
and  19U7  itself  is  important  as  the  year  in  which  George  B.  Dantzig 
developed  the  first,  and  so  far  the  best,  algorithm  for  the  solution  of 
the  general  linear  programming  problem:  the  simplex  method. 

Unfortunately  the  simplex  method,  as  it  is  used  on  computers  today, 
does  not  handle  every  problem  with  the  success  which  it  deserves.  At 
times,  with  no  forewarning,  computer  programs  which  implement  the  simplex 
method  will  produce  unreasonable  results  from  reasonable  data.  In  such 
cases  curses  are  muttered  against  a  mysterious  force  called  "round-off 
error",  and  a  battery  of  ad  hoc  panaceas  are  employed  to  attempt  a  cure. 
Sometimes,  in  hopes  of  forestalling  problems  from  the  outset,  a  variety 
of  heuristic  tests  and  purifications  may  be  built  into  a  computer  pro¬ 
gram  to  sniff  ,  at  suspected  indicators  of  round-off  error  and  attempt 
recoveries  during  the  computation.  But  all  of  the  medicines  employed 
have  the  aspect  of  folk  remedies.  They  are  largely  in  the  domain  of 
oral  tradition;  they  are  rarely  published,  and  their  efficacy  has  never 
been  proved. 

The  literature  concerning  the  effects  of  round-off  errors  upon 
the  simplex  method  is  very  small  (we  have  come  across  only  papers  [J], 
[8],  [10],  and  [  15 ]  listed  in  the  bibliography),  and  it  contains  no 
rigorous  results.  This  State  of  affairs  is  unfortunate  and  unnecessary. 
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June*  H.  Wilkinson  fir*t  published  the  tool*  for  investigating 
round-off  error  effect*  in  linear  algebraic  computations  in  1959,  [11], 
and  they  have  been  well-known  to  numerical  analysts  since  then.  Numer¬ 
ical  analysts,  however,  have  not  made  much  use  of  linear  programming. 

The  subject  is  Just  beginning  to  be  recognized  in  the  numerical  analyt¬ 
ical  world  as  a  piece  of  practical  machinery  (e.g.,  Rabinowitz  [9]). 

The  chief  users  of  linear  programming  up  to  now,  or.  the  other  hand, 
could  not  be  expected  to  be  familiar  with  the  literature  on  round-off 
errors  in  linear  algebraic  computations  (e.g.,  Chartres  and  Geuder  [2], 
and  Wilkinson  [11],  [12],  [13],  [14]).  Thus  a  proper  investigation  of 
round-off  errors  in  linear  programming  has  never  been  made.  This  thesis, 
we  hope,  will  serve  as  a  first  step  toward  changing  this. 

It  is  important  to  separate  the  abstract  expression  of  a  computa¬ 
tional  algorithm  from  the  realization  of  that  algorithm  on  a  computer. 

In  Chapter  2  we  present  the  simplex  method  in  an  abstract  setting,  and 
in  Chapter  3  we  describe  the  main  details  of  the  method' s  more  coranon 
implementations  on  computers.  We  show  in  Chapter  4  that  round-off  errors 
can,  indeed,  have  disasterous  effects  upon  the  conmon  simplex-method 
computer  programs.  This  is  not  due  to  any  flaw  in  the  simplex  method. 

It  is  rather  due  to  the  use  of  Gauss-Jordan  elimination  without  pivot 
selection  which  characterizes  the  standard  computer  implementations  of 
the  simplex  method.  In  Chapter  4  we  indicate  that  no  simple  a  priori 
inspection  of  the  data  will  yield  a  bound  on  the  effects  of  round-off 
errors  in  programs  which  employ  Gauss-Jordan  elimination  but  do  not 
permit  the  pivots  to  be  selected  arbitrarily. 


Chapters  6  and  8  present  alternative  computer  implementations  of 


the  simplex  method  for  which  pivot  selection  is  flexible,  and  these 
implementations  are  shown  in  Chapters  7  and  8  to  be  numerically  stable. 
That  is,  the  effects  of  round-off  errors  in  these  methods  can  be  bounded 
in  a  simple  a  priori  fashion  from  the  data. 

In  each  of  the  simplex -method  implementations  presented  in  this 
thesis  the  possibility  of  monitoring  round-off  errors  during  the  compu¬ 
tation  exists.  The  discussion  of  this  possibility  is  contained  in 
Chapters  5,  9*  end  10,  wherein  a  posteriori  error  bounds  and  their 
employment  in  simplex-method  computation  is  covered. 

In  rough  terms,  carrying  out  the  simplex  method  entails  the  repe¬ 
titive  solution  of  three  linear-equation  systems 

c 

Bu  =  b 
<  BTtt  -  y 
By  =»  A 

s 

w 

(The  notation  we  use  in  dealing  with  the  simplex  method  holds  closely 
to  that  used  by  Dantzig  in  [4].)  At  each  repetition  the  vectors  u, 
tt,  and  y  are  inspected,  the  basic  matrix  B  is  modified  as  a  result 
of  the  inspection  and  a  new  set  of  vectors  is  produced.  It  is  the  main 
effect  of  round-off  errors  that  only  approximate  vectors  u  +  6u, 
tt  +  6tt,  and  y  +  6y  are  produced  from  the  computations.  Assuming 
that  bounds  can  be  placed  on  the  norms  of  the  errors  ||6u||,  ||6tt||, 

and  ||6y||,  then  Chapter  10  discusses  the  use  of  these  bounds  during 


s implex -method  computation  to  make  allowances  for  round-off  effects . 
(Throughout  the  thesis  the  double  bars  |j  ||  denote  the  Lg  matrix 
and  vector  norms.  That  is. 


for  any  vector  v,  and 


wi  -  3s  & 


for  any  square  matrix  M  .) 

For  the  standard,  explicit-inverse  implementations  of  the  simplex 
method.  Chapter  5  covers  the  computation  of  error-norm  bounds .  The 
bounds  themselves  are  given  by  the  equation 

INI  <  ^-jj-E[[  ( I|e|| || v+fi v||+ ||f || || g|| ) 

appearing  in  mid-chapter,  where  fiv  represents  any  of  the  error  vectors, 
and  g  represents  any  of  the  right-hand  vectors.  A  bound  for  ||l||  is 
presented,  and  the  last  half  of  the  chapter  is  addressed  to  the  problem 
of  efficiently  computing  E,  or  at  least  an  upper  bound  for  ||e||  . 

The  other  two  implementations  presented  in  the  thesis  use  forms 
of  triangular  decompositions  for  matrices  as  substitutes  for  the  inverse 
of  the  basic  matrix.  The  standard  LU  implementation,  described  in 
Chapter  3*  dec  mposes  the  basic  matrix  into  triangular  factors  L  and 
U  .  The  modified  LU  implementation,  described  in  Chapter  6,  decomposes 
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the  basic  matrix  into  a  triangular  factor  L,  a  product  of  simple 
matrices  G ,  and  a  triangular  factor  U  .  It  is  shown  in  Chapters 
7  and  8  that,  for  both  implementations,  the  computed  solution  v  +  flv 
to  the  linear  system 


Bv  »  g 


satisfies 


(B+fi)(v+6v)  «  g  . 

Both  of  these  chapters  also  address  themselves  to  the  problem  of 
bounding  ||&||  .  For  the  modified  LU  implementation  the  results  are 
summarized  in  the  last  four  pages  of  Chapter  7*  For  the  standard  LU 
implementation  the  bound  is  given  in  Chapter  8  as 

llfcll  <  m[3rL+2m(m+J.)e'2]  lu^l  , 

where  the  quantities  and  are  defined  during  the  development 
of  this  bound,  and  m  is  of  the  order  B  . 

The  error  vector  8v  is  shown  in  Chapter  9  to  satisfy 

6v  =  -  B  ^6(v+6v)  , 


so  that 


||8v||  <  llB^llliellMvIl  . 
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Since  ||v+4v||  is  known  and  a  bound  for  ||£||  has  been  provided  by 

previous  chapters.  Chapter  9  addresses  itself  mainly  to  the  problem 

of  finding  an  upper  bound  for  easily  from  the  information  which 

ia  at  hand.  This  information  is,  roughly,  the  trace  and  determinant 
T 

of  BB  . 

T 

The  trace,  d,  of  BB  is  readily  computed,  but  only  an  approxi- 
T 

mation  to  det(BB  )  can  be  made.  The  analysis  of  Chapters  7  and  8 
show  that  the  computed  decompositions  of  B  used  in  the  standard  and 
the  modified  LU  implementations  actually  satisfy 

B  ^  fiB  =  decomposition 

because  of  round-off  er.or  effects,  where  a  bound  for  ||6b[|  is  given. 
Thus,  the  quantity 


p  =  det( [B+6B][B+6B]  ; 


is  at  hand.  It  is  shown  in  Chapter  9  that  the  best  upper  bound  on 
||B_I||  which  can  be  made  using  p  and  d  is 
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where  q  satisfies 


and  where 


K  >  ||  6  BBT+B6  BT+4  B6  BT||  . 

Finally,  in  Chapter  11  a  comparison  of  the  modified  and  standard 
LU  implementations  of  the  simplex  method  is  made  with  a  representative 
explicit -inverse  implementation  in  terms  of  speed  of  operation.  This 
is  done  by  counting  arithmetic  operations  needed  by  each  to  produce 
u,  tt,  and  y  and  to  accomodate  the  changes  in  B  which  follow. 

The  results  are  reviewed  on  the  last  page  of  the  chapter.  In  particular 
we  find  that  the  modified  LU  implementation,  which  we  have  shown  to  be 
not  subject  to  round-off  error  disasters,  is  competitive  in  speed  with 
standard  simplex -method  implementations. 


2.  Review  of  the  Simplex  Method  of  Linear  Programming 

To  begin  with,  we  will  briefly  cover  the  basic  concepts  and  termi¬ 
nology  of  linear  programming  and  of  the  simplex  method.  Everything  to 
be  mentioned  is  covered  in  greater  detail  in  Dantzig  [**)  and  in  Hadley 
[6]. 

A  linear  programming  problem  is,  by  definition,  one  which  can  be 
put  in  the  following  form: 

T 

(maximize  the  objective  function:  z  *  c  x 

subject  to  the  cor  straints  Ax  -  b  and  x  >  0  . 

A  is  an  m  X  n,  real  matrix  having  full  row  rank,  and  m  <  n  .  The 
vectors  b  and  c  are  real  and  have  dimensions  m  and  n  respectively. 

For  convenience  we  also  require  that  b  >  0,  but  this  is  net  essential 
tc  the  problem. 

In  the  standard  terminology  any  vector  x  which  satisfies  the 

constraints  is  a  feasible  solution,  and,  if  it  maximizes  the  objective 

function,  it  is  an  optimal  feasible  solution. 

Let  A  ,...,  A  be  ra  linearly  independent  columns  from  the 
1  m 

matrix  A  .  If  x  is  a  feasible  solution  having  x.  =  0  for  all 

J 

j  i  fv,  ,•••,  v_},  then  x  is  called  a  basic  feasible  solution,  and 

r  L  m  — -  . 

the  columns  A  ,...,  A  are  called  the  basic  vectors  (or  simply  the 
1  m 

basis)  associated  with  x  .  Any  m  X  m  matrix  composed  of  these  columns 
is  a  basic  matrix  corresponding  to  x  .  If  x^  =  0  for  some 
i  e  {v^,..«,  vm),  then  x  is  a  degenerate  basic  feasible  solution. 
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The  set  of  all  feasible  solutions  to  a  linear  programming  problem 
is  convex,  and  the  basic  feasible  solutions  comprise  precisely  the  set 
of  extreme  points.  The  existence  of  any  feasible  solutions  implies  the 
existence  of  at  least  one  basic  feasible  solution.  And,  if  the  objective 
function  is  bounded  from  above  on  the  set  of  feasible  solutions,  at 
least  one  of  the  basic  feasible  solutions  is  optimal,  (in  fact,  several 
of  the  basic  feasible  solutions  can  be  optimal.  A  given  linear  pro¬ 
gramming  problem  need  not  have  a  unique  solution. ) 

It  is  upon  bases  and  basic  feasible  solutions  that  Dantzig's  simplex 
method  is  grounded.  The  simplex  method  is  an  algorithm  which  searches 
from  basic  feasible  solution  to  basic  feasible  solution,  seeking  an 
optimal  solution.  The  search  is  designed  so  that  each  basic  feasible 
solution  encountered  has  a  greater  value  of  the  associated  objective 
function  than  its  predecessor.  The  search  proceeds  in  cycles.  Let  x 
be  a  nondegenerate  basic  feasible  solution.  The  simplex- method  cycle 
produces  one  of  the  following  from  x: 

I.  a  basic  feasible  solution  x'  from  which 


II.  an  indication  that  x  is  optimal; 

III.  an  indication  that  the  objective  function  is  unbounded  from 
above  on  the  set  of  feasible  solutions. 

Degenerate  basic  feasible  solutions  do  not  yield  quite  such  reasonable 
outcomes  from  the  simplex-method  cycle.  However,  there  are  standard 
ways  of  adjusting  a  given  linear  programming  problem  so  that  degeneracies 
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which  are  discovered  are  wiped  away,  which  fact  will  permit  us  the 
liberty  of  ignoring  degeneracies  from  this  point,  forward. 

To  each  basic  feasible  solution  there  is  a  corresponding  basis 
of  m  linearly  independent  columns  taken  from  the  n  columns  of  A  . 
And  outcome  (i)  of  the  simplex-metnod  cycle  guarantees  that  the  simple 
method  will  never  encounter  any  basic  feasible  solution  more  than  once 
Therefore,  starting  from  any  particular  basic  feasible  solution,  the 
simplex  method  cannot  be  carried  through  more  than 

n ' 

m.'  (n-m).' 


cycles  on  a  given  linear  programming  problem  before  terminating.  In 
practice  one  commonly  finds  that  the  number  of  cycles  taken  before 
termination  is  only  a  small  multiple  cf  m  (Dantzig  [U,  p.  l6o]). 

The  mechanism  for  producing  x'  from  x  by  the  simplex-method 
cycle  is  built  upon  basic  matrices  If  A  ,  A  is  the  basis 

V, 

1  m 

corresponding  to  x,  then  consider  the  basic  matrix 

B  *  (A  ,  - • •  ,  A  1  • 
v  -i  v 

1  m 

Let 

T 

u  ■=  [x  ,  . . ,  X  ] 

V-,  v_ 


be  the  vector  of  the  nonzero  components  of  x  ■  Then 


KVrr  <  r'W  i 


Bu  =  b 


To  carry  out  a  simplex-method  cycle,  one  begins  with  B  and  proceeds 
as  follows: 

1.  Solve  Bu  *  b  • 

T 

2.  Let  V  *  [c  c  ]  ,  the  vector  composed  from  c  by 

1  m 

selecting  those  components  corresponding  to  the  basis. 

T 

3-  Solve  B  tt  =  y  . 

U.  For  each  column  Aj  of  A  not  in  the  basis  compute 


cj  =  "  Aj  '  V 


and  choose  any  ind£x  s  for  which  c  <  0 

s 


5-  Solve  By  =  A  • 
s 

6.  Let  r  be  an  index  such  that 


Ur  min  _^i 

yr  ‘  y>o  y4 


7-  Drop  the  vector  A  from  the  basis  and  add  the  vector  A  . 

V  s 

r 

The  resulting  set  of  m  vectors  is  a  new  basis,  whose 
associated  basic  feasible  solution  is  the  vector  x7  .  The 
associated  basis  matrix  B7  can  be  constructed  with  whatever 
ordering  of  columns  is  most  convenient. 


The  cycle  cannot  b;  carried  past  step  3  if  all  Cj  are  nonnegative, 
in  which  case  x  is  optimal.  And  the  cycle  cannot  be  carried  out  past 
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V 


ctep  5  if  no  component  of  the  y  vector  is  positive,  in  which  case 
the  objective  function  is  unbounded  on  the  set  of  feasible  solutions. 
The  change  in  the  value  of  the  objective  function 

•  T  /  /  \ 

z  -  z  =  c  (x  -  x; 


is  given  by 


r  /  T.  \ 
—  (c  -n  A  )  = 
y  s  s' 
r 


Thus  it  makes  some  sense  to  choose  s  such  that 

c„  =  min  c.  , 
s  j  J 

making  c  as  largely  negative  as  possible,  in  hopes  of  making  as  great 
s 

a  change  in  the  objective  function  as  possible. 

To  solve  a  given  linear  programming  problem  by  the  simplex  method, 
a  basic  feasible  solution  must  first  be  found.  If  none  exists,  the 
problem  is  infeasible.  The  task  of  finding  the  initial  basic  feasible 
solution  can  be  formulated  as  an  auxiliary  linear  programming  problem, 
one  for  which  a  starting  basis  is  known.  Hence,  the  simplex  method  is 
generally  carried  out  as  a  two-phase  process,  each  phase  of  which 
consists  of  repeated  applications  of  the  simplex-method  cycle.  The  main 
part  of  this  thesis  will  treat  with  the  computational  details  of  carrying 
out  the  simplex-method  cycle  on  an  automatic  digital  computer  and  with 
the  round-off  errors  incurred  during  the  computation. 
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).  Standard  Computer  Implementations  oi1  the  Simplex  Method 


The  major  price  to  be  paid  in  carrying  out  the  simplex-method 


cycle  is  incurred  by  solving  the  linear  systems 


T 

B  TT  =  v 


By  -  A 


If  these  systems  were  presented  in  isolatior  we  would  expect  to  spend 

3 

an  order  of  m  arithmetic  operations  on  a  computing  machine  to  get 
their  solutions.  In  the  context  of  the  simplex  method,  however,  this 
much  work  is  not  usually  required.  The  basic  matrix  B  used  in  one 
execution  of  the  simplex-method  cycle  is  constructed  from  all  but  one 
of  the  columns  which  formed  the  basic  matrix  used  in  the  preceding 
execution  of  the  simplex-method  cycle.  The  vector  b  remains  constant 
throughout,  and  the  vector  Y  changes  in  only  one  component  at  each 
cycle.  These  observations  provide  t'  0  opportunity  for  significant 
savings  in  computational  effort. 

For  example,  in  representative  computer  implementations  of  the 
simplex  method  the  matrix  of  order  m  +  1 


1  -Y 

D  =  0  B 


c  constructed-  »' j  1 : 1  < •  e 


-1 


1 

'V 

V  ‘ 

0 

i 

it  can  be  verified  that 


**  — 

-  --  D* 

u 


0 

b 


1 


1 

1 

1 

1 


and  that 


-  - 

l  . 

n  A  -c 

3  S 

'  s 

-  I)-1 

-c 

s 

y 

M  a* 

y 

A 

Let  us  number  by  zero  the  initial  components  of  the  vector s,  and 
tne  initial  rows  and  columns  cf  the  matrices  in  what  follows  directly 
below •  This  will  unify  the  indices  used  in  this  chapter  with  those 
used  throughout  the  rest  of  the  tnesic. 

Suppose  an  execution  of  the  simplex -met hod  cycle  ends  in  outcome 
( I ) .  If  we  let 


D'  = 


i 

•T* 

0 

B' 

14 


] 


i 


by  replacing  the  r  column  of  D  with 


then  it  can  be  veril’ied  that  (D' )  ^  is  obtainable  from  D-^  by  the 
application  of  the  Gauss-Jordan  elimination  which  reduces  the  vector 


r 


where  e^  is  the  rth  unit  vector,  then 
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(D'  )"]  =  TD'1 


Similarly  we  can  update  the  vector  u,  rather  than  solve  a  system  of 
equations  1'or  u  : 


In  this  manner  the  simplex-method  cycle  car.  be  carried  out  on  a  computer 

2 

with  only  an  order  of  m  arithmetic  operations. 

If  D  1  is  explicitly  computed  in  an  implementation  of  the  simplex 
method,  is  kept  iri  computer  storage  as  an  m  x  m  matrix,  and  is  updated 
by  Gauss-Jordan  elimination  at  each  cycle,  the  implementation  is  said 
to  use  the  explicit  inverse-  Alternatively,  D-1  can  be  expressed  as 
a  product  of  transformations  which  reduces  D  to  the  identity  matrix. 

At  the  end  of  each  execution  of  the  simplex-method  cycle,  a  Gauss- 
Jordan  transformation  can  simply  be  added  to  the  product  so  that 

D_1  =  , . .  -r^1)  . 

Implementations  of  the  simplex  method  built  along  these  lines  are  said 

to  use  the  product  form  of  the  inverse.  They  are  used  when  the  data 

matrix  A  is  sparse  and  has  its  nonzero  entries  in  a  regular  pattern 

(i  ) 

which  permits  the  use  of  transformations  T  '  having  a  particularly 
simple  form. 

The  basic  computational  process  in  implementations  of  either  type 
is  the  Gauss- Jordan  elimination  applied  to  a  vector,  either  a  column 
of  D_1,  the  vector 
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» ■  '»ljJ 


‘*n» . 


I 

I 

I 

B 

n 

n 

n 

n 

n 


or,  in  case  the  product  form  of  the  inverse  is  used,  to  the  vector 


Hence,  we  will  want  to  study  this  process  as  it  is  carried  out  in 
floating-point  arithmetic  on  contemporary  computing  machines. 


0 

I] 

D 

0 

0 

0 

C 

B 

I 

[ 


i 


.  Hound-off  Krror  AnuLysis  u‘  Guuss-Jordan  Elimination 


We  will  denote  the  computed  value  of  an  arithmetic  expression  l 
obtained  by  floating-point  operations  on  a  computer  by 


l't(Z) 


This  value  is  a  function  not  only  of  l  but  also  of  the  computer 
employed  and  the  algorithm  used  for  the  evaluation.  Following  Wilkinson 
[lj],  we  know  that  on  reasonable  computers,  sc  long  as  overflow  or 
underflow  do  not  occur, 


fl(a*b)  = 
l'/(a-b)  * 
ff(axb)  = 


a(l+H1)  +  b(l+Hg) 
a(UT15)  -  b(l+H4) 
a  x  b  x  (1*1^) 


ft  (a/bi  «  (a/b)  x  (1*116) 


for  any  two  machine  numbers  a  and  b,  where  the  quantities  |f|^|  are 
bounded  by  some  fixed,  small,  positive  number  In  fact  we  will  assume 
that  our  computations  are  carried  out  on  a  machine  offering  a  normal- 
precision  floating-point  arithmetic,  for  which 


hi  I  <  •• 


and  an  extended-precision  floating-point  arithmetic,  for  which 


5  £,  ■ 


I 


where  c is  a  few  orders  of  magnitude  smaller  than  •  This  is 
usually  accomplished  by  carrying  out  floating-point  arithmetic  in 
base  t  arithmetic  (t  «*  2,  8,  10,  l6  are  common),  carrying  a 
significant  figures  in  normal  arithmetic  and  Og  significant  figures 
in  extended  arithmetic  (usually  with  Og  >  2o^)  •  If  that  is  the  case, 
one  may  take 


e 


1 


T 


1-0. 

1 


and 


A  normal-precision  machine  number  is  a  number  which  may  be 
represented  exactly  on  the  computer  to  the  precision  of  normal  arith¬ 
metic.  A  similar  definition  holds  for  extended-precision  machine 
numbers. 

Suppose  we  transform  the  vector  (w^,...,  w^)  into  the  vector 
(v^,...,  Vj)  by  applying  Gauss- Jordan  elimination,  obtaining  the 
multipliers  from  the  vector  (p^,...,  p^),  with  p^  used  as  the 
pivot  (l<p<l)  .  The  p's,  v's,  andw's  are  all  to  be  normal-precision 
machine  numbers. 

Thus, 


vp  *  > 

vi  =  f|(wi"  Piwp/V  *  f°r  1  ^  P  * 
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Or 


ve  “  (wB/pp)(14/n)  » 

v i  -  wi(l+cpi)  -  (piWp/Pp)(l+Tl)(l-Hii)(l-»pi)>  for  i  i  3  . 


This  means  that,  letting  flwQ  =  T\w0  , 

P  P 


'b  ■  ?“  (W  • 


And,  for  i  i  e  > 


'i  ‘  l"i-  5;  V 


p 


Or,  letting 


wit*>i- 


we  have 


Ti  ■  (W  •  r  • 

P 


Thus,  if 


P 


[I 


±  (p-e»))eto)  ] 
*P 
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it  hue  been  estubli shed  Unit 


v  -  u(Pw)  =  p(w+6w) 


That  is,  the  round-off  errors  have  the  effect 
vector  v  be  an  exact  transform  via  tne  given 
of  a  perturbation  on  the  vector  w  . 

How  large  can  the  pertumv  i  on  flw  be? 
Unless  restrictions  arc  placed  upon  the  pivot 
6w  cannot  be  bounded.  We  have,  for  i  *  0  : 


of  making  the  computed 
Gauss-Jordan  elimination 


The 


answer  is  disturbing, 
the  components  of 


IP.  I 

flWjl  <  IwJItpJ  +  jj-j  |wp|  j^.+p^.p.+TU.+Tlp.-UU^.I  , 


or 


K1  i  |uilei  '"e112"!*  3ei  +  ‘i>  ■ 


And  so,  if  | p  |  is  made  small  enough  relative  to  the  magnitude  of 

P 

some  one  of  the  other  components  of  the  vector  p,  the  right-hand 
side  of  the  latter  inequality  can  be  made  arbitrarily  large  for  the 
corresponding  value  of  i  . 

To  see  a  trivial  example  of  what  this  means,  suppose  that  we 
have  a  computing  machine  at  our  disposal  which  performs  decimal  arith¬ 
metic  with  tnree  figures  of  precision  in  noxmal  floating-point  calcu¬ 
lations.  The  result  of  each  basic  arithmetic  operation,  however,  is  to  be 
computed  to  six  figures  of  accuracy  before  being  rounded  to  the  nearest 
three  significant  figures-  This  means  that 


2: 


t, 

1 1 


£1  B  I10’2 


w  =  (1,1)  ,  p  =  (10"  ,1)  ,  and  0=1. 


v1  =  fi(l/lO-1+)  -  lo\  and 

Vp  -  fi(l-lxi/l0"4)  =  -io"u  . 


That  is, 


10"  0  1 


It  4 

v2  -10  -10  1  0 


So  fiw2  ■  -w2  .  A  perturbation  has  been  made  in  one  of  the  components 
of  w  which  is  as  large  as  the  component  itself.  If  there  was  any 
useful  information  in  the  vector  w  initially,  it  has  been  lost. 

It  is  easy  to  concoct  examples  wherein  relatively  small  pivots 
arise  during  execution  of  the  simplex  method.  E.g.,  take  the 


(1  1  2] 


1  0  Y  [y/2 

“  ,  b  - 

1°  1  y  llJ 


I 


where  V  is 
columns  of 


and 


Since 

proceeds- 


a  smull  positive  .lumber  less  than  one.  The  firs,  -wo 
L  will  serve  as  an  Initial  ousts,  giving 


r 

Z 

"1  ♦  T/2" 

U1 

- 

Y/2 

.U2- 

1 

_  - 

ttTA,  -  c,  =  Y  -  1»  whlch  is 
;>  ■> 

In  the  notation  of  Chapter  2, 


less  than  zero,  the  cycle 
s  =  5,  and  we  must  compute 


Booh  components  of  y  are  positive 


und 


I 


3 

] 

I 

I 

I 


u 

y 


1 

1 


1 

2 


Sc  —  -  min, 


basic  vector 


or  r  =  1 


in  the  notation  of  Chapter  2. 


Hence,  the 


must  be  exchanged  for  the  basic  vector 


Thus  all  necessary  updating  is  to  be  carried  out  by  the  Gauss-Jordau 
elimination  composed  from  the  vector 


using  y^  =  Y  as  the  pivot,  which  is  small  relative  to  the  other  tvo 
components. 

The  implication  of  all  of  this  for  standard  implementations  of  the 
simplex  method  is  that,  since  the  relative  sizes  of  the  pivots  which  arise 
during  computation  cannot  be  easily  predicted,  no  bound,  simply  expressed 
in  terms  of  the  data,  can  be  made  a  priori  on  the  perturbations  which 
will  be  caused  by  round-off  errors. 


2h 


M4U> 


5.  Error  Bound;;  l  or  LI  near- cyst  cm  Liolutior.r 
Found  by  Explicit  Inverse- 

Though  it  i s  not  (ossible  a  priori  to  bound  round-off  perturbations 
during  application  oi’  the  simplex  methcc  computer  implementations 

such  as  those  just  discussed,  we  cun  compute  round -oi’f  error  bounds  in 
an  a  posteriori  fashion.  We  will  do  that  here  in  connection  with 
explicit-inverse  implementations.  Product- foim  implementations  can  be 
treated  in  a  similar  fashion. 

Consider  the  following  problem:  if  F  is  a  given,  nonsingular 
matrix  of  order  t,  if  g  is  an  i  vector,  and  if  the  matrix  Q  is 
accepted  as  an  approximate  inverse  to  P,  how  much  error  is  maoe  in 
taking 


f/(Qg) 


as  the  solution  to 


Pv  =  g  ? 

Let  v  =  P  and  v  +  6v  =  ff(Qg)  •  Each  component  of  v  +  6v 
is  formed  as  an  inner  product: 


\  +  !vi  “ f,<£  <V'o)  • 


The  calculation  of  an  inner  product  is  so  basic  to  computational  linear 
algebra  that,  in  order  to  be  as  accurate  as  possible,  it  is  a  frequent 
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practice  to  form  the  products  and  partial  sums  in  extended-precision 
arithmetic.  Only  the  full  sum  is  reduced  to  normal  precision  (Wilkinson 
(13,  P*  23]).  That  is, 


t<°>  -  0 

t(J)  =  ♦  qjjgjfUSjKl+Tlj) 


for  J  ■  1  , . . . ,  I  , 


where  l®jl»  IPjl,  and  iTljl  are  all  bounded  by  c g  .  Only 
reduced  to  normal  precision,  giving: 

vi  +  #vi  -  • 


where  (pj  <  •  Therefore 


'i  +  8vi  ■  (l*Pi ) (l^j ) U-tlj  ^  <1<0k> 


or 


l 

\  ■  ,vi  *  - 


where 


is 


2  6 


V 


Thus,  in  matrix  form, 


v  +  6v  =  (Qrt>)g  , 


where 


(i+tj 


L±1 

t 


[exp( 


)-l)} 


The  bound 


el  +  ^lf£i)  '~f^2  iexP(  i~-  )-lJ 


it. 


>  |(l+p  )(l*S  )(l+T|.)  f[  (1+cr  )-l 
1  J  il 


k=j*l 


will  be  established  in  Chapter  7 
Let  the  residual  matrix 


E  *  I  -  QP 

be  defined.  We  will  demand  that  Q  be  such  a  good  approximation  to 
the  inverse  of  P  that  ||e||  <  1  .  Then  the  inverse  of  I  -  E  will 
exist,  and 


II(I-E)’1||  < 


SOT  • 
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Therefore 


v  +  6v  =  Qg  +  'tg 
=  QPv  +  <fg 

-  QP(v+fiv)  -  QPfiv  +  <fg  , 


which  gives 


5v  =  (I-E)  ~l<i>g-E(v+6v) ]  , 


or 


IM  <  x.jlElf  (llE||||v+6v||+||t||||g||)  . 


The  cost  of  computing  the  residual  matrix  E  is,  in  ordinary 

situations,  far  too  great  to  pay  for  simply  putting  a  bound  on  ||6v||  . 

3 

An  order  of  f  arithmetic  operations  is  needed  to  form  E  .  That  is 

essentially  the  same  amount  of  work  as  is  required  to  invert  the  matrix 

P  .  There  are  far  more  lucrative  ways  to  expend  one's  energy.  The 

often-used  technique  of  iterative  refinement  (Wilkinson  [13,  p.  130], 

Forsythe  and  Moler  (5,  Chapter  13]),  for  example,  can  correct  v  +  6v 

to  a  point  at  which  ||6v||  w  and  can  be  carried  out  in  an  order  of 
2 

l  operations.  But,  in  the  context  of  the  simplex  method,  the  com¬ 
putation  of  this  error  bound  can  become  quite  economical. 

Suppose  P7  is  obtained  from  P  by  replacement  of  the  rth  column 


Pe 


(r) 


2-S 


with  the  vector  w  : 


P'  -  P  -  (Pcr;-  w)e  r 


Let  the  vector  t  be  given  by 


t.  =  f  /  ( Qw )  , 


and  compute  Q* ,  an  approximation  to  (?)  >  by  apj lying  to  Q  the 

Gauss- Jordan  elimination  which  reduces  t  to  the  r  unit  vector: 


Q'  =  fi{[l  -  i-  (t-c^r^)c^r^  jG.} 


(This  imitates  the  situation  in  typical  explicit-inverse  implementations 
of  the  simplex  method. )  Apply  the  results  cf  Chapter  L  columnwise  to 
Q  to  get: 


Q'  =  [I  -  ~  (t-e^r '  )e^  ] I QtY ]  , 


where  for  all  i  and  j  : 


If  .  .  |  < 

*  ij  - 


—  q  .  E.  ti  1  -  r 

t  rj  1 

r 


O  el+  |-U|  |,  |  (2,1+5.^5)  if  i^ 


For  convenience  let 


T  =  [I  - 


f  <t-e''r))e(r)  J 

r 


i 


Then 


E'  =  I  -  qV 

=  I  -  T(Q+y)p' 

■  I  -  TQP'  -  TfP'  . 

But 

I  -  T»'  -  I  -  [I  -  f  (t.e(r))e(-')Tl4[P.(Pe(r).  v)e(r^] 

r 

-  1  -  («■  -  f  (t-e(r))e(r>V«(Pe(r)-  „)e(r)T 

r 

+  F  (t-e^)e^  Q(Pe^-  w)e^r^] 
r 

*  E  -  f  *£“  (t-e^)e^  Q(Pe^-  w)e^T 

r 

-  (t-e^Je^  QP-Q(pe^r^-  w)e^r^T]  * 
r 

We  make  the  substitution 


t  «  fl(Qw)  =  (Q*l)w  , 


for  which 


<  {e1+  (l+e1) 


t-A+2 


fexp( 


111 

1-c, 


)-!])  k 


ij 
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and  the  substitution 


QP  “  I  -  E 


to  obtain 


I  -  TQP'  =  E  -  (  j-  (t-e^')(e^  e^)e^ 

v 

r 

-  f  (t-o(r))e(r)TEe(!')«!(r)’i 
r 

.  i.  (t-o(r))(e(r)Tt)e(r;T 


.  i  (t-e(r))e(r)‘ 
r 


+  f  (t-e(r))e(r)TE  -  e(r)e(r)' 


+  Ee(r)e<r)Tt  te(r)T} 


t  nve(r)T  -  f  (t-e(r))e(r)Tnve(r)T 

r 


A  couple  of  simplifications  are  possible: 


f  (t-e(r))(e(r)Ta(r))a(r)T  »  f  (t-e(r))e(r)I  , 
r  r 

since  e^  e^r^  -  1  , 


and 


t 

r 


(t-f»)(e(r) 


t)- 


(0J 


since 


t  . 
r 


(t-c(r))e(r) 


Therefore 


T 

I  =  TOP'  =  S  -  {  7-  (t-e^r^)e^r^  E(l-e 
r 

+  Ee^e^  }  +  mwe^ 

=  [I  -  (t-e^)e^  ]E(l-e 

r 

(r)T 

+  mveKrj  . 


(r)e(r)T} 

(r)e(r )T) 


So 


E'  =  T[E(l-e^r^e^r' )  •  Owe^  -  YP' ]  . 

And  therefore, 

He' II  <  ||t||  (||E||+||n|| IIwII+IIyIHIp'  || )  . 

If  upper  bounds  for  ||e||,  and  ||y||  are  at  hand,  an  upper  bound 

for  He'II  is  thereby  available.  And  only  an  upper  bound  for  the 
residual  matrix  norm  is  needed  to  bound  the  error  J|fi v||  . 
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Returning  to  the  notation  and  concept:;  of  Chupt<rs  2  1  i»  i.* 

clear  that  the  above  results  can  bo  employed  dui-irip.  t  ,e  <•><.  w  !  \ 

the  simplex  method  cycle  to  bound  the  computational  error  irici  -rod  ir 
finding  u,  y,  tt,  z,  and  the  •  In  a  la‘ei  chapter;  in  con:v  :ti 
with  another  implementation  ol'  the  simplex  method;  we  will  cv.^ri  ier  li  e- 
employment  of  such  error  bounds  to  provide  automatic  tolerance  corral. 
The  techniques  presented  there  can  easily  be  imitated  with  the  explicit, 
inverse  implementation  of  the  simplex  method. 


I 

I 

1 

1 

i 

7 

-* 

I 

1 

I 

1 


***** 


r 


I 

I 
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6.  The  Mocill'led  LU  Implementation  of  the  Simplex  Method 

Returning  tc  the  notation  of  Chapter  2,  w c  let  B  be  the  m  x  m 
nonsingular  basic  matrix  and  concern  our. solver.  w",h  salving  the  linear 
systems 


...» 


m 

Bu  =  b  ,  B*rr  =  Y  .  By  =  . 


-»♦ 

I 

-  • 

•  • 

! 

•  • 

T 

-  » 


Suppose  we  effect  this  by  decomposing  B  into  the  product  of  a  lower 
and  an  upper  triangular  matrix.  Such  a  decomposition  is  always  possible 
provided  that,  at  the  least,  row  permutations  of  B  are  permitted. 

Column  permutations  could  be  added,  but,  in  tne  context  of  the  simplex 
method  in  which  we  work,  it  will  not  be  convenient  to  permute  the  columns 
of  each  basic  matrix  we  encounter  arbtrarily. 

Thus  we  have 


IlB  =  LU 


*? 
•  • 

j 

1 

I 

I 

I 

I 


for  some  permutation  matrix  n,  upper  triangular  matrix  U  and  lower 
triangular  matrix  L,  where  we  further  demand  that  all  diagonal  elements 
of  L  be  +1  . 

A  given  linear  system 


Bv  =  q 


which  we  must  solve  is  handled  by  solving  the  triangular  systems 


Lt  ■  II  q 


Uv  ■  t  . 


The  round-off  error  analysis  of  the  decomposition  and  solution  of 
the  triangular  systems  is  given  in  Wilkinson  [13].  The  best  error  bounds 
are  obtained  if  n  is  chosen  so  that  all  elements  of  L  are  bounded  in 
magnitude  by  unity.  The  strategy  of  constructing  II  during  the  com¬ 
putation  of  the  decomposition  has  come  to  be  called  "partial  pivoting" 
or  "row  pivoting". 

At  the  beginning  of  the  simplex  method  cycle  B  has  the  form 


B  =  [A  ,  •  •  • ,  A  ]  • 


At  the  end  of  the  cycle  we  may  be  required  to  construct  a  new  basic 
matrix  B^^  from  the  columns 


A  , « •  * ,  A  »  A  ,  •  •  • ,  A  ,  A 
V1  Vr1-1  vr1+l  vm  1 


And  this  is  to  be  done  in  such  a  way  that  the  transition  from  the 
decomposition  of  B  to  one  for  B^  is  particularly  easy  and  com¬ 
putationally  stable. 

Consider,  then,  the  basic  matrix 


B  =  [A  ,  •  •  • ,  A  }  A  , . . . ,  A  ,  A  ]  ; 
V1  vrL-l  \+l  vm  S1 


W* 


I 

I 

I 

I 

I 

I 

I 

T 

1 

7 

-  » 

-  • 

I 

—  • 

'! 


.1 

w9 


/  *  \ 

i.e.,  let  B  be  obtained  from  B  by  dropping  the  r  rojurm, 

shifting  all  subsequent  columns  one  position  to  the  left,  und  appending 

A  on  the  right.  Therefore,  columnwise: 
s. 


iT'W1)  =  [L_1nA  L_IT1A 


rri 


L  V  >  • « • ,  L  ^TIA  ,  L  ^TlAr  ] 
vr^+l  vm  U1 


[u, ,  u  , ,  u  U  ,  L’"T1A  ] 

1  rr1  ri x  m  si 


=  Hll>  , 


,(1) 


where  the  U.  are  the  columns  of  U  .  ir  '  is  a  matrix  of  the  form 

l 


1 

I 

I 

1 

1 


that  is,  an  upper  Hessenberg  matrix  with  zeros  below  the  diagonal  in 

was 

produced  as  a  byproduct  during  the  computation  of  the  vector  y  .  Hence 
can  be  constructed  without  any  computational  effort. 


the  first 


r  -1 
1 


columns.  We  point  out  that  the  vector  L 


~hA 


3  6 


I 


I  I  m m  ******  «***- 


^  can  be  reduced  to  an  upper  triangular  matrix  by  u r>int; 

Gaussian  eliminations  to  zero  out  the  subdiagonal  elements  in  columns 

through  m  .  The  partial  pivoting  technique  is  carried  out  by  making 
a  decision  at  each  elimination  step  whether  or  not  to  interchange  two 
adjacent  rows.  Thus,  is  gotten  from  by  applying  a  sequence 

of  simple  transformatins  on  the  left: 


uW.AV1} ...  , 

m-1  m-1  r,  r. 


1  ‘1 


where  each  has  the  form 


i  i+1 


1 

• 

*  1 

1 

* 

-g(1) 

Ki 

1 

1 

• 

• 

'  1 

. 

.  i 


and  each  is  either  the  identity  matrix  or  the  identity  with  the 

ith  and  i+lSt  rows  exchanged,  the  choice  being  made  so  that  |gp^|  <  1  • 
Thus 


IIB 


(1) 


nir1)r^1) 

i  i 


’1  ...  n(1]r(1^V1> 

m-1  m-1 
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2 

I 

I 


mrnmmmmmmmm 


Thisi  is  n  decomposition  or  B 


U) 

(!) 


ultublc-  lor  our  purpose:;.  We  could 


(id, just  me  order  ol'  row.  in  B  arid  multiply  out 


in(  0r( 1  ) 

ri  1 1 


-i 


m-l  m-J 


to  produce  a  decomposition 


noyn  .  ,(!!,/:) 


which  copies  the  decomposition  used  1'or  B  .  But  this  turns  out  to 

require  more  work  1'or  the  execution  ol'  the  simplex  method  cycle  than 

is  required  if  the  decomposition  is  left  in  the  product  form  given  above. 

If  another  execution  of  the  simplex  method  cycle  is  taken,  the 

resulting  basic  matrix  B^  is  to  be  formed  from  ^  in  the  same 

( 1 ) 

manner  in  which  B  was  formed  from  B  .  The  decomposition  for 

(2) 

B  which  we  would  use  ? n  the  subsequent  cycle  would  therefore  be 


nB<2) 


ri  U 


m-l  m-l 


X  [n(2)r(2)']  ...  n<2>r(2!'V2) 

r2  rp  m-l  m-l 


And  in  general,  if  is  the  k^h  basic  matrix  which  we  construct, 

the  general  decomposition  which  will  concern  us  will  be 


•T> 


AWdtwiuw"  . .  ■  " 


iib^  =  un^M1) 

ri  ri 


m-l  m-1 


-1 


(nU)r(k) 

r,  ri 
k  k 


-1 


n(k)r(k) 

m-l  tn-1 


Each  given  linear  system 


B 


q 


is  solved  by  computing  the  solution  to 


Lt  =  Ilq  , 


carrying  out  the  transformations 


w 


rrk>"rk)'  -  tr<lMlJ 

rR  rk  m-l  m-l 


VL^]t 


> 


and  solving 


In  following  chapters  we  shall  show  that  the  round-off  errors 
committed  in  the  above  computations  can  be  bounded  a  priori,  which  was 
not  the  case  in  the  standard  implementations  of  the  simplex  method. 

And  we  shall  show  that  the  computational  effort  required  for  the 
modified  LU  implementation  compares  favorably  with  the  effort  required 
to  carry  out  a  representative  explicit-inverse  implementation. 
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7 .  Hound-off  Minor  Analysis  of  the  Modified 
LI1  Imp  lenient  at  i  in 

Wilkinson  (12,  Vj  |  has  analysed  the  comj.ur.ut  ional  error  associated 
with  solving  the  linear  system 

13  v  -  j 

via  the  LU  decomposition  of  B  as  fellows: 

I.  It  is  assumed,  without  loss  of  generality,  that  the  rows  of 
B  have  been  prearranged  so  -.hat  the  partial  pivoting  strategy 
will  net  alter  their  erdei i ng • 

II  The  matrices  L  and  U  which  are  computed  from  B  satisfy 


B  -  6B  =  LU 

III.  In  solving  the  triangular  systems  composed  from  L  and  U, 
vectors  w  and  t  are  produced  which  satisfy 


I 

I 

I 

I 

I 

1 

I 


(L*5L)w  =•  q  , 

( U+6  U )  t  -  w 

And  t  in  the  computed  approximation  to  v  . 

IV.  Bounds  are  given  for  |j6B]|,  i|5L||,  and  ]|6Uj|  •  Thus  the 

vector  t  has  been  shown  to  satisfy 


•mmkuw «. wart^aw*™***"**"' 


(B+C)L  =  ( D+fi B+L6 ' 1+6 LU+fi L6 U ) t  =  q  , 


with  a  known  bound  for  II fill  . 

We  wish  to  establish  a  similar  result  for  the  decomposition 


b^  =  L[n^r^  ...  n^V1]  ]  ...  [n^kVk^  ...  n^kMk|  ]u'k^  , 

1  r,  r,  m-1  m-1  r,  r  m-1  m-1 


which,  for  notational  convenience,  we  will  also  write 


B(k)  .  J^uyk) 


The  following  error  analysis  will  be  shown  to  hoid: 

I.  It  is  assumed  that  the  arrangement  of  rows  in  the  first  basic 
matrix  B  is  such  that  no  row  rearrangement  is  needed  to 
produce 


(B+5B)  »  LIJ 


computationally  by  the  partial  pivoting  strategy. 

II.  At  the  simplex-method  cycle  the  matrices  L  and  l/k^ 

(k) 

and  the  product  G  satisfy 

B^k^  +  6B^  =  . 

4i 


I 

I 

I 


III.  In  solving 


I 

I 

I 

I 

I 

I 

I 


B'^v  =  q 


via  the  decomposition  above,  vectors  f,  w,  and  t  are 
produced  which  satisfy 


(L+6L)w  =  q  , 
(G^+flG^V  =  vi  , 


und 


] 

I 

I 

I 

1 


(u^+6l/kV  =  f 


And  t  is  the  computed  approximation  to  v  . 

( lr  }  !  Ir  ^ 

IV.  Bounds  for  the  norms  of  the  matrices  6B'  ,  6L,  flU'  , 


fiG 


00 


are  given.  Thus  the  vector  t  will  satisfy 


(B(k)+£(k)jt  =  q  , 


where 


=  fiB*^  +  LG^kVd^k^  +  L6G^l/k^ 

+  6LG^l/k  ^  +  L6G^6l/k^ 


+  fiL&G^W*^ 


and 


42 


(k) 

Hence,  a  bound  for  ||£'  /||  is  known. 

Chartres  and  Geuder  [2]  give  the  following  useful  lemma: 

If 

ICj  <  e  <  1 

and 

0  <  p,  <  v  <  a  , 

then 

H  v 

|  n  a+Ci)  f[  (l+Ci)*1_ll  <  ^[exp(  ^  )-l]  . 

i=l  i=u+l 

They  also  implicitly  give  the  corollary: 

If 

tp  < 

i  n  u-tyj  n  <  m  , 

j=i  j  j=P+i  j 

and 

V>  v 

|  1]  (l+c,  )  n  U+Ci)"1-!!  <  N  , 

i=l  1  i=p+l 

then 

p  k  pi  v 

i  n  d+cpj  n  (i+cpj_i  n  (1+Ci)  n  (no-1-!!  <  m + n + mn 

j=l  J  j=p+l  J  i=l  i=p+l 


1 


I 

I 

I 

I 

T 


t 

\ 

-to 


t 

-to 


7 


] 


I 

I 

I 

I 

I 


This  is  useful  if  the  bounds  on  the  jcpjl  differ  from  these  on  the 


ICil 


We  proceed  to  justify  the  claims  made  under  II,  III,  and  IV  just 
preceding. 

Dropping  the  superscript  (k)  momentarily  for  convenience,  we 
review  the  Chartres-Geuder  demonstration  that 


(U+6u)t  =  f  , 


-1 


where  t  is  the  approximation  to  U  f  computed  by  back-substitution 


The  components  of  t  are  produced  in  order  from  t  down  to  t.. 

m  1 


by  the  computation 


m 


with  partial  sums  accumulated,  as  usual,  in  extended  precision.  That 


is, 


Pi  - 


»2  *  V1+V  - 


m-i 


m 


m-i 


rm-i+l  =  fi  H  ( 1+\ )  •  t  li ,|t/(1+>A|.i)(1+CTi.i)  n  (L+\) 

v=l  x=i+l  v=l-i+l 


14 


! 

I 


7. 


Finally, 


where  |l|  |  ,  luj  ^!»  and  |<7 ^  .  |  ure  all  bounded  by  . 

ti  =  (l+cp)(l+p)  , 

i,i 


where  cp  and  p  are  the  errors  due  to  the  reduction  of  p  ...  to 
^  y  *n-i+i 

normal  precision  and  to  the  division  by  i  ,  respectively.  Thus, 

1 ,  i 

m-i 

.t^d+cpj-^i+p)"1  ]]  (i+n  y1 

v=l 


.  'l1 

*  ^  <lV  . 


or 


f 


i 


m 


m 


L  ui  1*1  +  L  *ui  t*i 

i= i  ifl  1  /=i  1 


> 


where 


[(l+cp)"1(l+p ) 


m-i 


> 


and 


i-i 


!ui,«  -  n  (nrl-ilui,i 


for  l  >  i  •  Appealing  to  the  lemma  and  corollary  just  given,  we  may 
set  the  following  bounds: 


where 


2c  -t? 

1 6ui ^  i  |  <  {(m-Oeg  +  - 2  ( l+(m-i Jeg]}  ^  ui , i 


lftUi,jl  <  ^(j-i+2)£2^Ui,jl  f°r  j>i 


'2  m+2 


lf  ({m 

+2  ^ exp  l  ~ 


(m+2)e0N 


•2  J 


’-1}  • 


But  using  the  inequality 


.  max  | .  | 

i,j  Ka1  ’ 


which  holds  for  all  matrices  D  (Wilkinson  [12]),  we  may  write 


2E  -e2 

lltull  <  »((«*!)£'  [l+(ra-l)Eg]J  “j  |utiJl  • 


In  the  same  manner  it  is  shown  that 


(L+6L)w  =  q  , 


,<“>  =  f<ir(k,)n.lk.)w(tt'-1) 


m-1  m-1 


r  =  „(a>)  , 


where 


k 

u>  =  Z  (m~ri )  • 

i=l 


Thus  we  wish  to  make  an  analysis  of 

d  =  fl(ma) 


for  any  vector  a,  where 


and  where  II  is  either  the  identity  matrix  or  the  identity 
Vth  and  v+l3^  rows  interchanged. 

If  II  is  the  identity,  then 


d  =  fl(lTIo) 


is  equivalent,  compenentwise,  to  the  following: 


a.  for  i  1  v  +  1 

\  ' 

=  fVrgav]  +  tav+1O-gav(0+T+pT)] 


where  |o|,  |p|,  and  |t|  are  all  bounded  by  £ ^  .  Thus, 

<P  =  -g(p+t+PT)  , 


with  the 


if  we  let 


we  have  in  matrix  form  d  =  (r+fif)a,  where 


or 


d  ■  (rn+«m)a  . 


So  again 


d  -  (r+6r)na  . 


Thus,  for 


f  -  fi{r£k]n£k|  ...  r^n^w) 

m-i  m-i  r.  r. 


we  may  write 


f  "  [rikl+ari1Cllnikl 

m-l  m-i  m-i 


[r^+firS.1)  ini1^ 


or 


n(k>[r(k}+«r(kh- 

m-l  m-l  m-l 


Therefore 


where 


0<k>  -  ntl)r<l>‘1 ...  n^A)'1  , 

r,  r,  m-1  m-1  ’ 


1  ‘1 


and 


eoc<) .  na)cr(ir^rarl] 

ri  ri  ri 


Recall  at  this  point  that 


n<k>(r<k>'Vk>'1).G(k)  . 

m-1  m-i  m-i 


np>n<J) 


for  all  possible  i  and  J  .  If  we  interpose  such  products  at  appro- 

(k) 

priate  places  among  the  factors  in  $G  >  we  can  arrange  to  express 
(k) 

6G  in  the  equivalent  form 

6G(k)  -  p^W^W1)]  ...  ...  , 

r^  r^  m-1  m-i  r^  m-1 


where 


00 

m- 


i 


>'+ 


II 


•nd 


°nkl  +  6nikl  “  rikl  1  +  6r^kj  1  > 

m-1  m-1  m-1  m-1 

°mk2  +  "  nikl[rik2  L+6rmk2  > 

m-z  m-2  m-l  m-2  m-2  m-1  ' 


n* ^  +  fin^  .  p(^[r^1^  +6r^  . 

rl  rl  rl  rl 


(k) 

This  puts  the  product  6G'  into  a  form  to  which  we  can  suitably 
apply  the  following  lemma: 

If  H^, ...»  11^  and  G^, . . . ,  Gp  are  square  matrices,  then 


«L 


a  •  • 


<H1-°1)02 


•  •  • 


G 

P 


+  H^Hg-G^Gj  ...  Gp 

+ . . .+  ...  Hp_^( Hp“Gp)  • 
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wwwwwu.  I  .'HA.milAWIIWMIiatWMi 


According  to  this, 


60(k)  «  ...  nf kj i 

r,  r, m-i 


1  ‘1 


+  ((n^^W^)  ...  (o^+6n^kh«r/kh} 

r.  r,  m-c  m-e  m-i 


1  *1 


Each  has  diagonal  elements  equal  to  +1  and  only  one  non¬ 

zero  off-diagonal  element,  gp\  whose  magnitude  is  bounded  by  one. 

This  element  appears  somewhere  in  the  subdiagonal  portion  of  the  matrix: 


1 

\ 

!(J)\ 

h  K 

1 

The  corresponding  6n^J'  has  only  two  nonzero  entries: 


= 

li)  J 

U) 

5 

L  "i 

• 

_ 1 _ 

V<J)  4J) 
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Hence 


f 


. 


IUo(lt)l|  <  |||#0(l,)|||  <  I|p(k)||  ||  E  E  (iS^Mnf0!) 

1  J 


K 


But 


!|p(k)ll  -  (In*15  n^ll  »  i  , 


and 


<  m  , 


and  there  are 


E  (m-ri) 

i  * 


i-1 


(k) 

terms  in  the  expression  for  fiG'  .  Therefore 


^e,+e,  k 


/  .  \  hE.te  ft 

||60  k  II  <  m  -jij-i  £  <-rt)  • 


Finally  we  wish  to  establish  the  relationship  between  the  basic 
matrix  and  the  factors  L,  G^),  and  of  its  numerically 

computed  decomposition.  We  proceed  inductively. 


fQ 


■■■ . ...  .  „ 
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fl 

0 

D 

0 

D 

D 

0 

C 

D 

0 


:  1  iimiiiiaiiini' 


mmmv  Wfl  •  vm&*rtnm 


If  B^  is  the  basic  matrix  initially  decomposed,  then  matrices 


L  and  IT  are  computed  such  that 


B^  +  fiB^  *  Ll/°^  • 


For  the  moment  we  neglect  the  superscript  (0)  .  The  decomposition 
is  produced  by  letting 


=  f'<V  E.  'i*V  f°r  ‘SJ 

v=l 


'id  ■  fl[(V  Ex  for  i>J> 


where  J  is  assigned  in  order  the  integers  from  1  to  m,  and  for 
each  value  of  j,  i  is  assigned  in  order  the  integers  from  1  to 
m  .  The  inner  products 


i-1  J-l 

V  t ,  u  and  Y  l,  u  . 

v*l  iv  V=1  iv  vJ 


are,  as  usual,  to  be  computed  in  extended-precision  arithmetic.  Thus, 
for  each  pair  of  values  for  i  and  j,  if  we  let 


H  =  min{i,j}  , 


60 


then  we  set 


pl 


> 


and 


P  *  p  (1+T)  )  -  i,  u  ,(l+a  )(1+C  ) 
v+1  v  v  iv  vj  v  *v 

for  v  *  1»  2,.m,  |i  -  1  , 

where  |l^|,  |a  |»  and  I  Cv  I  are  bounded  by  e2  •  Finally, 
if  i  <  j  we  set 


=  p(1+cp) 

P 


f 


and  if  i  >  j  we  set 


‘i 3  ~  VUJJ)(1^)(l+p)  ’ 

where  |<p|  <  represents  the  error  committed  in  reducing  to 
normal-precision  accuracy,  and  jp|  <  represents  the  error  committed 
while  performing  the  division.  We  may  also  write,  for  the  sake  of 


Uy  =  P^(l-hp)(l+p )  > 


Di. 


symmetry. 


V 


iJWVJWWfW  !  ►vwfc-'m  I 


1 

1 

I 

I 


il’  p  =  0  in  this  case.  Thus,,  in  general, 


VVj  “  u+<p)(i+p)pu  , 


where  we  have  used  the  fact  that  *  1  .  We  find  also  that 


1 

1 


M.=l 

n  <i+\> 

V=1 


I 

3 

3 

3 

3 


]_  |JL  —  1 

■tww  n  (»v  • 

<--l  V=K+ 1 


Therefore 


U-l 

b« -  vn  d+\> 


K=1 


'iAj 


d«K)(ucx) 


n  ^v1 

V=1 


3 

3 

1 

I 

1 

I 


That  is, 


b 


ij 


£ 

K=1 


liKUKi 


> 


6? 


■ 

■ 


i 


where 


r~ 

4bu  ■  n  d^)'1-1) 


1  V— 1 


2e i+e? 

|Sbij|  s  |litiui»j|l(“'1),2 +  ^Tf  d+d-D^n 


And,  since 


we  have 


+  Y.  I  liKUKj^K+2^e2 

K=  1 


|6B||  <  l.  ™  l«bul  , 


2  +  2 

||6B||  <  m2  ju,  Jf.(ra+l)e2  +  — - — “2  t1+(m+1)e2^  * 

^  (l-e^ 


Now  suppose  that,  after  K  exchanges,  we  have 


)(k)  +  6B^  =  LG^kVk)  • 


That  is,  columnwise 


[B(k),...,  Bik)l  +  iB^k): 


-  LG(k)[U^k),-..,  U^]  . 
i  m 


Let  B^k+1^  be  formed  as  described  in  the  preceding  chapter.  That 

A,}  fk}  at 

is,  drop  column  B'  '  from  Bv  ,  shift  the  r.+l  through  the 
th  K 

m  column  left  one  place,  and  append  the  column  A  .  The  deccm- 

*k  (k) 

position  is  made  to  reflect  this  change  by  modifying  U  .  The 
(k) 

column  IT  '  is  dropped,  the  subsequent  columns  are  shifted  one 
k  / »  , •> v 

position  to  the  left,  and  the  column  IT*  ,  which  is  the  computed 
approximation  to 


g^Va  , 

sk 


is  appended- 


r(k+l) 


From  previous  discussions  we  know  that  Hv  '  satisfies 

m 


(L-tA^)(G^kWk))H(k+1^  =  A  , 


where 


||a^||  <  { (m+2)tg  +  [l+(m-l)cg]}m  , 
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. n 

. 


'ffammmrmVttkt.UtiKWI-smimMi*  r.  MW, 


and 


/ 1_\  4e,+ef  k 

""  llin"^ri(n'r‘ 


Therefore 


LqOc)^!)  s  .  [^^(^(jOc^kykJjndcfl) 


m 


A.  +  6a 


Sk  Bk 


So  that  columnwise 


,(k) 


b'«1(  •«.....  b<b),  v 

+  (44k)’ . Bvr  «B<k>,  .A.  ] 

»B(kH>  +  e(kH> 


m  '  s. 


“(k)[uik)’->  u^r  »<k),  H'k+l)) 


=  LG^^H^k+i^ 


At  this  point  the  transformations 


r(kfl)n(k+l)  (k+1)  (kfl) 

1  m-l  a-1  *••  rr,  nr. 


,  •» 

■  a  & 
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V 


are  applied  to  to  produce  .  The  product  of  transfor¬ 


mations 


G(k)n(k+l)r(k+l)_1  _  n(k+l)p(k+l)* 
rk  rk  m-1 


.(k+1) 


is  taken  as  G 

Columnwise  we  have 


(k+D  =  R(k+i) 


(k+i)  _  H(k+i) 


IT  T'  =  H 

V1 


V1 


and 


u(k+l)  .  n(!r(kn)n(k+i)]H(kn), 


ri  ri 

k  k 


i/1?1’  -  «( [rik:1)nik:1>  ... 

m-1  m-1  m-1  r.  r,  m-1  * 

k  k 

u'k+1>  .  fiur^^ni11:1’  ...  4k+1>n<k+1>)H<k+1>5 

m  m-l  m-l  r,  r.  m  J 


Hence,  from  previous  discussions, 


(k+1)  _  „(k+l) 


H,  =  U. 


H 


(k+1)  _  y(k+l) 


V1 


V1 


oo 


for  each  j  =  1,  . . 


m  by 


HU+1)  .  ^k+D^k+l)  4  y(k+l) 
J  m  j 


whe  re 


v<ktl>  .....  tC"})  .  0 

1  V1 


and 


y(  k+i)  =  g^k+^^k+l) 


v  v 


for  v  =  r^, ...»  m 


Thus 


H(k+1)  =  *(k+1)u(^l)  +  y(k+l) 


m 


where  Y^k+1^  is  the  matrix  with  columns  Y^k41^ 


.  Therefore 


B^k+1^  +  ©(k+1)  = 


=  LG^k^[<t^k+1^U^k+1^+Y^k41^  ] 


m 


Or,  if  G^k+1^  =  G^kVk+1^ 

m  ’ 


,(*!)  +  e(k+l)  _  LQ(k)?(k+l)  _ 
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« 

i 

r 


That  is, 


B(*n)  +  JB(ku)  .  m(k*i)0(k-M)  ( 

and  the  induction  is  complete. 

To  bound  ||6B^k+^||,  we  note  that 


>(k+l)  =  [ftB^,...,  fiB^,  flB^,...,  ®Ag  ]  , 

k  k  k 


6A  =  -  [LE^k^^k^G^k^-»A^k^H^k^]H^k+1^  , 
sk  m 


y(k+l)  a  [y(k+l)  <<>  y(k+l)  y(k+l)  y(k+l) -i 

1  1  >"  ’  >  r^_i  *  >'  •  •  >  ra  J 

.  [0,...,  0,  S»(k+l)u(^)] 


r,  r, 
k  k 


m  m 


II® A  ||  <  (llLlI  ||H(k)||+|^(k)||  j|G(k)||+||A^k)||  ||E(k)||)||Hi5k+l)| 


2  4el+el  £ 


<  m  -5-7 —  j  (m-r.)+{(m+2)e'  +  77-  [l+(m-l)e']} 
E1  i=l  *  i"el 


4e,+e?  k  "1  \ 

t1+T^rI  (»«1)3jiiH‘kl)i 

i  1-1  j 


which  we  denote  by  .  So,  if 


0(k)  >  max  |8b(k)|  ( 
~  1  ^  J  J-J 


C {/S(k),  P(k*l}}  >  max  |o.  . 

—  1>J 


Furthermore 


max  |  (K+l)|  <  max  y^kU^^Di 
ijJ  ~  J  J  <3 


<  raax  i|6$(k+l)||  maX  |ju(k+l)| 

J  J  r  r 


2  kel*El  ,  ,  max  i  (k+l)i 

<m  ___  (mH-rk,  i;j  lujj  'I  , 


which  we  denote  by  Therefore,  we  may  conclude  that 


max  |8b(kU)!  <  max  |9(k*l),  +  m  max  |  (k*l),  ( 
AjJ  1J  -LJ  1J 


and  that 


||6B^k+1^||  <  ||  0  ^k+1^||  +  ||LG^||  ||Y ^k+1^|| 
<  m[max{/?^,  P'k+1h  +  mR^k+1^] 


eh  "V 


To  recapitulate:  If  we  are  asked  to  compute  the  vector  v 

satisfying 


B 


OO 


v 


■  q 


via  the  modified  LU  decomposition 

+  6B^  =  L(j(k)u(k)  f 

we  actually  calculate  the  vector  t,  an  approximation  to  v,  which 
satisfies 


(gOO+gMjt  =  q  , 

where 

=  6B^  +  LG^k^6U^k^  +  L6G^k^U^k^ 

+  &LG^k^U^k^  +  L6G^fil/k^ 

+  fiLG^k^fiU^k^  +  6L6G^l/k) 

+  6L6G^kWk)  . 


Now 


l|6L||  <  m{(m+2)( 


[1_ 

-e. 


[l+(m-l)eg] 3  , 
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V 


/u\  Jn-l+e,  Jl 
l|6G(k)||  <  m-ji— ±  t  » 


•1  i=l 


2  « 

ll«U(k)||  <  mt(m-l)e'  +  ~~2  tl+(m-l)e']l  “  |  .'k)| 


And  in  order  to  bound  || 6 Bv  '||,  we  let 


s(0)  = « iu<°>n(m+i)e'  +  flilii, 


2  [l+(ra+l)cg]) 


e<k+1>.ma*tsHp<k+1^+n*(k+l! 


for  k  =  0,  1,  2,...  .  This  makes 


||#B(k)||  <  m8(k) 


for  all  k  . 


Thus  it  is  possible  to  compute  a  bound  on  j|£v  'll  at  every 
execution  of  the  simplex-method  cycle.  This  will  be  the  foundation 
for  the  a  posteriori  bounds  developed  in  Chapter  9*  Moreover,  an 
a  priori  bound  on  ||&^||  for  all  k,  expressed  only  in  terms  of 
the  given  data,  could  be  given  at  this  point.  The  bound  would  be 
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extremely  large  and  of  no  practical  value,  but  its  existence  argues  a 
numerical  stability  for  the  modified  !U  implementation  of  the  simplex 
method  which  was  not  present  in  the  explicit-inverse  implementation. 


This  a  priori  bound  would  be  grounded  upon  the  following  considerations, 

(k)  (k) 

round-off  effects  in  the  computation  of  u^  '  and  h ^  being  neglected 

for  the  moment: 


k  < 


n-' 

m.’  (n-m)  I 


> 


m  -  r.  <  m  -  1  for  all  i  , 
1  — 


max  iu(0)i  <  pm-l  max  ,a 
i<j  ij  -  i,j  ij 


(Wilkinson  [13,  p.  97])#  and 


max 

i<j 


<  m 


max 

i#j 


h 


(k) 

ij 


for  all  k  =  1,  2, . . . 


(Wilkinson  [14,  p.  218]),  where 


max  | h(k)  | 

i#  j  ij 


r  max 
<  maxi  .  , 
-  ifJ 


m 


and 


H^||  =  ||(G^k'+-^k^)_1(L-»A^k^)'1A 
m  £ 


k-1 


<  m  ?ax  la. .1  (approximately). 

-  i,j  1  ij1 


Finally  we  must  note  that  one  of  the  systems  of  equations  to  be 

(k)T 

solved  in  the  simplex-method  cycle  uses  B  '  as  the  coefficient 
matrix.  It  can  be  shown  that,  if  t  is  the  solution  to 


*<*>  v  .  , 


computed  via  the  modified  LU  decomposition  for  ,  then  t  satisfies 


(B(k)  +S'(k))t  =  q  , 


where  any  bound  on  ||ft  ||  is  also  a  suitable  bound  for  ||^ 


.-.I 
i  ■+ 


8.  The  Standard  LU  Implementation 


of  the  Simplex  Method 


The  round-off  error  analysis  due  to  Chartres  and  Geuder  L 2 j  which 
holds  for  the  conventional  method  of  solving  a  linear  system 


Bv  =  q 


via  the  LU  decomposition  of  B  was  covered,  in  passing,  throughout 
the  discussion  in  the  preceding  chapter.  It  is  worth  a  quick  review 
for  its  own  sake,  since  we  wish  to  describe  an  implementation  of  the 
simplex  method  in  which  this  error  analysis  is  valid. 

The  decomposition  of  B  is  produced  by  letting 


u.  .  =  ff(b  .  -  Y  t.  u  .)  for  i  <  j 

ij  lJ  v=i  lv  VJ 


'i.i =  •  t  'i„Vi )/«,,)  for  1  >  ^  • 


iJ  V=1  lv  VJ  JJ 


From  this  we  find  that  L  and  U  satisfy 


B  +  6B  *  LU  , 


where,  letting  y,  =  min{i,j]  , 


2  +  2 

l^ijl  -  +  ~ — ~o  U+k-lJeg]) 


(1-6  J 


K=1 


A  vector  w  satisfying 


(L+4L)w  =  q 


is  produced  by  letting 


i-1 

wi =  f,(tli  -  'uV 


for  i  =  1,  m,  where 


Ifi^iil  <  (1-1)63  +  TT^  [i+(i-i)«2l  > 


and 


Ifi^ijl  <  Uijl(J+2)«2  for  1  >  J  • 


Finally  a  vector  t,  approximating  v  and  satisfying 


(U+*U)t  =  w  , 


7b 


V 


is  produced  by  letting 


t. 

1 


r':(ui  -  £  v^i1 


for  i  -  m,  m  -  1 , . . . ,  1,  where 


|6i  |  <  f  (m-i)Eg  +  — — -2  [ l+(m-i )tg ])  |i  u|  , 

(1-6]) 


and 


I fiui j !  5  (J"i+^)£2luij  I  for  J  >  1  • 


Hence  the  vector  t  satisfies 


(B+e)t  =  q  , 


where 


6  =  6B  +  6LU  +  L6U  +  6L6U  , 


or 


(61.  u  ,+t.  6u  .+61.  6u  .) 

IK  KJ  IK  KJ  IK  K  j  ' 
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So 


ieiji  s  ivvjII(i*'1)£2  +  7IT7-)2  u+ta-u^n 

+  l*V"nJl  +  lV“Wl  +  '‘VVl1 

+  ^  UlKUKj|((j+K+6)t2+(K+2)(3-K*2)Ej2)  , 


where,  if  we  let 


\  = 


2t.+t? 

(m-l)eg  +  - 0  [l+(m-l)£^]  , 


(i -«1r 


we  may  sst  the  bounds 


I'Wj1  -  < 


M-J 


Xluyl 

for 

i  <  j 

X|UJj' 

for 

i  =  J 

(■HlUijU.jU' 

for 

i  >  i 

(m+1) | u±J | eg 

for 

i  <  J 

Xlujj 

for 

i  =  J 

^Vjj1 

for 

i  >  j 
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and 


XOm-lJIUjjItg  for  i  <  J  , 


"W  -  {  x2,ujj' 


for  i  =  j  ,  and 


xOn+lJli^  -j j  |  for  i  >  J  • 


This  gives  the  result 


M.-1 

I j  I  <  I  lLKaKj  I  ^  (.i+/<+6)  +  (/<+2)  (.1-K+2)e2]e2 


k=  1 


f  [2X+(m+l)e2+(m+l)Xeg]|uij|  for  i<j 


+  \  [3X+\2]|u^|  for  i  =  j 

[2X+(m+l)e2+(m+l)XE2]|i.Ju^|  for  i  >  j 


Let  us  set 


=  e2(l+2e2) 


and 


=  X  +  (m+1 )eg  +  (m+l)e2X  +  X2  • 


As  estimates,  then,  we  will  have 


~  ' 

^  e2  ^  e2 


and 
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and 


8o 


can  be  computed  within  each  execution  of  the  simplex-method  cycle.  This 

quantity,  being  only  a  constant  multiple  of  |u^|  *  ^as  better 

behavior  than  the  computed  bound  on  llfc^il  developed  in  the  preceding 

chapter.  For,  initially,  the  modified  LU  implementation  and  the  standard 

LU  implementation  are  the  same.  That  is,  has  the  bound  which 

we  have  just  given  above.  On  subsequent  executions  of  the  simplex- 

(k) 

method  cycle,  however,  terms  depending  upon  the  transformations  G  ' 

(k) 

enter  the  bound  for  ||&'  ;||  .  These  terms  bring  in  the  quantities 

£  (m-r^)  and  /3^  , 

i=l 

which  grow  roughly  in  a  linear  fashion  with  k  and  cause  the  bound  for 
||&^)||  gradually  to  grow  larger  than  the  bound  given  above  for  ||g.||  . 

Finally,  we  should  note  as  in  the  preceding  chapter  that  if  we 
solve  a  linear  system 

T 

B  v  =  q 

via  the  LU  decomposition  computed  for  B,  then  we  actually  calculate 
a  vector  t  which  satisfies 

(BT+?)t  =  q  , 

where  any  bound  on  ||&||  is  also  a  suitable  bound  for  j|?||  . 

The  error  analysis  given  above  can  be  used  in  the  context  of  the 

Cj 


simplex  method  as  follows. 


Let 


B  -=  [A 

V1 


be  the  current  basic  matrix,  and  let  L  and  U  be  the  matrices  of 
its  computed  decomposition.  We  assume  that  the  order  of  the  rows  in 
B  is  such  that  the  partial  pivoting  strategy  produces  no  row  reordering. 
If  the  simplex  method  cycle  is  carried  out  to  produce  the  new 

basis 


r-1 


r+1 


m 


As> 


then  let  the  new  basis  ms.tr ix  B  be  constructed  by  dropping  the  r 

column  of  B  and  introducing  A  in  its  place.  In  this  case  the 

s 

decomposition  matrices  L'  and  U*  which  would  be  computed  from  B; 
using  the  partial  pivoting  strategy  differ  from  L  and  U  only  in 
the  shaded  areas  indicated  below: 


r 


l 


Thus,  only  in  the  shaded  areas  must  the  elements  l[  .  and  u/ . 

ij  iJ 

be  computed,  and  only  for  the  last  m  -  3  +  1  rows  do  row  interchanges 
have  to  be  considered-  The  saving  of  work  implied  by  this  is  considered 


in  detail  in  Chapter  11.  But  in  rough  terms,  half  the  work  needed  to 

compute  L  and  U'  from  scratch  is  needed  on  the  average  to  update 

the  LU  decomposition  of  B  to  the  L'U'  decomposition  of  B'  -  Even 

3 

so,  this  amounts  to  an  order  of  m  operations  for  the  average  execution 

of  the  simplex  method  cycle-  The  conventional  implementations  of  the 

simplex  method  as  well  as  the  modified  LU  implementation  described  pre- 

2 

viously  require  only  an  order  of  m  operations- 


We  have  seen  that  the  round-off  errors  committed  in  computing  the 
solution  v  to  the  linear  system 

Bv  =  g 

permit  only  an  approximate  solution  v  +  6v  to  be  produced  which 
satisfies 

(B+£)(v+4v)  =  g 

for  some  matrix  t  whose  norm  we  can  bound.  Thus 

flv  =  -  B-16(v+6v)  , 

and  so 

INI  <  Kb"1)!  ||e||  ||v+M|  .• 

Since  ||v+fiv||  can  be  computed  and  an  upper  bound  on  ||e||  is 
available,  we  can  bound  ||6vj|  by  finding  an  upper  bound  for 


B  +  «B  =  VW  , 

where  det(v)  and  det(w)  are  particularly  easy  to  compute,  and  a 
bound  for  ||6b||  is  available,  then  we  can  compute  an  approximation  to 
the  product  of  the  eigenvalues  as  well: 

P  =  [det(B+6B)]2 
m 

=  n  x. ([B+fiB][B+fiB]T) 
i=l 

m 

=  ]]  X.  (BBT+6BBT+B4BT+6B«BT)  . 
i=l 

But  BBT  and  6BBT  +  +  6B$BT  are  both  symmetric  matrices.  Hence, 


it  follows  that 


Xi(BBT+6BBT+BfiBT+4B«BT)  =  ^  +  6\ 


for  each  i,  where 


1 6X  L  |  <  ||6BBT+B4BT+6B4BT|I  . 

(See  Wilkinson  [l4,  Chapter  2,  §44].)  Suppose,  therefore,  that 

K  >  ||«BBT+B6BT+6BSBT||  . 


Then 


IfixJ  <  K 


for  all  i  .  And 


m 

p  =  n  , 

i=l  1  1 


m  in 

d-  £  -  £ 


6Xj 


i=l 


i=l 


The  best  lower  bound  for  \±  =  Xmin(BBT)  which  we  could  hope  to  obtain 
using  this  information  would  be  the  component  of  the  solution  to 
the  following  programming  problem: 


I 

tl'Whmw  <■  -  ■  « 


r 


minimize 


subject  to 


f(V  —  ’  V  ’  qi 


m 


81(q1»-**»  qm)  *  Y  qi  -  d  -  mK  <  0 


g2(q1i-*- »  q^)  ®  Y  ^  "  d  +  -  0 

i=l 

m 

6j(qi>*  •  •  »  \)  =  Y  1°e(qi)  -  log(p) 


=  0 


Using  the  method  of  Lagrange  multipliers,  we  are  presented  with  four 


linear  systems  to  solve: 

1. 

r 

m 

Y  lo8(q, )  -  i°g(p)  =  ° 

i=l  ± 

wA-° 

-  U1/q2  =  0 

•  • 

•  • 

•  • 

*  •  ‘‘A  * 0 

which  has  no  solution; 


[ 

| 

i 


m 

Y  losK )  -  l°g(p)  =  0 

i=l  1 

m 

q.  -  d  -  mK  ~  0 

l 


1  '  ^l^ql  "  ^2  =  0 

-  ^1/q2  -  =  0 


•  “A  -  t*8  " 0  • 


from  which  we  first  infer  that 


q2  =  ...  =  ^  -  Q  > 


giving  the  equations 


m-l 


P 


(m-l )Q  =  d  +  mK  , 


which  imply  that  q^  must  satisfy 


/  m-l 

ql  =  d+mK-q. 


Nm-1 


3. 


log(qi)  -  log(p)  =  0 


m 


i  1  -  M>1/q1  -  = 


0 


0 


-  jij/qa  -  =  0 

•  • 

•  • 

-  “A.  ■  ■  0  ’ 


Ki 


WQVfMmmivrti- 


whose  solution,  as  with  the  preceding  system,  has  a  first  component 
which  must  satisfy 


/  m-1  vm-1 

ql  =  d-mK-q.  ^  ; 


]T  los(<3i )  -  log(p)  =  0 
i=l 


y  q,  -  d  -  mK  =  0 
i=l  1 


Y  q,  -  d  +  mK  =  0 
i=l  1 

1  -  -  M-2  -  =  0 

"  P>]/ ~  ^2  ”  ^5  =  ^ 

•  • 

•  • 

•  # 

-  '  ^2  -  ^3  =  0  ' 


from  which  we  first  infer  that 


q2  "  qm  ~  Q  ' 


giving  the  equations 


r  nm_1- 

•---  p 


^  q^  +  (m-l)Q  =  d  +  mK 
ql  +  (m-l)q  *  d  -  mK  , 


which  have  no  solution. 

Thus  we  use  the  lower  bound 


\±  >  rnaxtt^,  q_}  , 


where  q.,  q  satisfy 

f  .  "Nm-l 
m-1  I 

^  =  p^d-hnK-qf  J 

If  F  is  a  real  matrix  of  order  m  with  positive  eigenvalues, 
the  fixed  point  f  of  the  equation 

f  ■  ’"'1 

is  often  used  as  a  lower  bound  for  Xmin(F)  (Bodewig  [l,  page  69]). 
The  fixed  point  is  found  iteratively  with  f  =  0  inserted  initially 
into  the  right-hand  side.  With  this  background,  it  is  easily  shown 
that 
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f: 


|  I 

!  i 

1 

i 

i 

i 

i 


i 


n 


or 


,(ih)  „  (in) 


-  '»+ 


Hence 


q  =  lim  q^  >  lim  q^1^  =  q,  . 

,  —  ,  +  + 
l  -*  oo  i  -<  » 

Therefore,  if  Z  >  ||fc|j,  we  can  bound  Jjfi v||  by 
IMI  Z||v+6v||  . 


:i 


] 

i 

l 


I 

I 

I 

I 

I 

1 

.1 

'? 

-4 

7 


n 

n 


:1 

il 

!l 

H 

3 

3 

1 
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10.  Employment  of  Error  Bounds  in  the 
Simplex-Method  Cycle 

Each  study  which  we  have  made  on  an  implementation  of  the  simplex 
method  has  led  us  to  conclude  that,  in  solving  the  systems 


T 

B  TT  =  v 
(  Bu  =  b 

.By  -  As 


computationally,  round-off  errors  permit  us  to  produce  only  the  approxi¬ 
mate  vectors 


tt  +  fiir,  u  +  6u,  y  +  6y  • 

From  these  studies  we  have  also  determined  how  bounds  on  the  errors 

||6tt||  <  u> 

<  M  <  ^ 

Jlfiyll  <  cp 

can  be  computed. 

Furthermore,  in  computing  the  objective  function 

T 

z  =  tt  b  , 
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wi:  can  produce  only 


z  +  6z  = 


,  T 

(n+6n  j  c 


-  n^'b  +  6n'Lb 


z  +  6tt'  Ii  , 


so  that,  letting 


5  =  > 


|flz|  »  | fin"1!1 1  <  ||4tt||  ||b||  <  % 


Likewise,  in  computing 


C.  =  T!  A.  -  0. 

J  J  J 


l’or  j  /  tv.,...,  v  },  1  <  j  <  n,  we  can  only  produce  c.  +  6c.  , 

m  J  J 


where 


|6c  |  =  |6ttTA  |  <  £  . 

J  J  J 


for  £  =  oj||A  .|| 

«J  J 


These  bounds  may  be  used  to  give  warnings  that  round-off  effects 
have  become  harmful.  As  an  example,  at  each  execution  of  the  simplex- 
method  cycle  the  test  that  B  is,  indeed,  a  basic  matrix  is  that 
u  >  0  .  Conceivably  round-off  errors  committed  during  the  previous 
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execution  of  the  simplex-method  cycle  cun  have  caused  th  ?  wron  •  cojuun 
to  be  moved  to  or  from  the  basis.  As  a  check  on  u.,  then,  \/c  note  that, 
for  any  i, 

u,  >  0  if  u.  +  6u,  >  T]  , 

l  l  l  1 

and 

u.  <  0  if  u.  +  6u.  <-  V  . 

i  ii 

If  neither  of  the  above  hold,  then  we  shouldn't  draw  conclusions  about 
the  sign  of  u^  .  Thus,  if  for  any  i  it  is  discovered  that 

ui  +  6ui  <  T1  > 

then  the  computation  in  the  previous  execution  of  the  simplex-method 
cycle  bears  a  second  look. 

Similarly  an  improper  exchange  of  columns  in  the  previous  execution 
of  the  simplex-method  cycle  may  have  produced  a  basic  matrix,  but  one 
for  which  the  associated  objective  function  value  shows  a  decrease 
rather  than  the  increase  which  it  should.  Let  z  +  Jz  and  £  be 
respectively  the  previously  computed  objective  function  value  and 
associated  error  bound.  Let  z  +  6z  and  5  corresponding 

values  for  the  current  execution.  Then 

z<z  if  z  +  6z+£<z  +  6z  +  5, 


*&>«***«! 


and 


z  <  z  if  z  +  6z  +  §<z  +  6z+|. 


If  neither  holds,  the  situation  is  indeterminate.  Thus,  if 


z  +  6z  +  5<z+6z+|  f 


the  computation  of  the  preceding  simplex-method  cycle  could  stand  review. 
Likewise  the  index  s  for  which 


c  =  minfc 
s 


has  definitely  been  found  if 


c  +6c  +C  <c.+6c. 
s  s  ®s  —  j  j 


-  c 


j 


for  all  j  j.  [v,,...,  v  )  and  1  <  j  <  n  •  And  c  is  definitely 
x  m  *“  s 

negative  if 

c  +  ic  +  C  <0. 
s  s  *s 

The  execution  of  the  simplex-method  cycle  can  proceed  if  a  value  of 

s  is  found  such  that  only  this  last  inequality  holds.  Otherwise  each 

c.  is  nonnegative  or  has  indeterminate  sign,  and  as  good  a  check  as 
J 

possible  should  be  made  on  the  optimality  of  the  current  basic  feasible 


solution. 


9o 


1 

I 

I 

I 

** 

I 

I 

-l 

i 


k. 


.  i 


IT  there  ii>  no  value  of  j  for  whici 

y.  f  £y.  >  cp  , 

J  K> 

then  each  y.  is  negative  or  has  indeterminate  sign,  and  as  r,oo<!  a 
J 

check  as  possible  should  be  made  on  the  sign  of  each  component,  of  y 
to  establish  whether  or  not  the  objective  function  is  unbounded  from 
above.  If  for  each  j  the  sign  of  y.  can  be  determined;  i-c., 
ei ther 


y,  +  6y.  >  cp 

J  J 


or 


y,  +  6y ;  <  -  tp  , 

J  J 


then,  if  r  satisfies 


3 

1 

J 

1 

1 


u  +Su  +n  i  .+6u  .-7] 

£  r  1  <  j  j  1 

y +3yr-'P  -  .  i+6y,+cp 

i  -  J  J 


for  all  j  such  that  y.  +  5y.  >  cp,  r  is  definitely  an  index  such 

J  0 

that 


u 

r 


min 


u . 


is  any  index  for  which 


V01 
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Thus,  tests  involving  Tj,  U),  and  cp  can  be  made  throughout 
e:  jcution  of  the  simplex-method  cycle  to  determine  whether  t!vj  compu¬ 
tations  which  have  just  been  performed  need  to  be  investigated-  The 
needed  investigations  must  be  carried  out  by  producing  better  approxi¬ 
mations  to  the  vectors  u,  tt,  and  y  corresponding  either  to  the 
current  basis  or  to  the  proceeding  one,  and  computing  smaller  error 
bounds  T),  ur,  and  cp  in  association.  Triis  is  most,  easily  accom¬ 
plished  by  iteratively  refining  the  three  vectors  (Wilkinson  [15, 
pp.  121-127],  Forsythe  and  Moler  [5,  Chapters  13  and  22],  Moler  [?])• 

Under  the  usual  circumstances;  i.e.,  if  e0  is  sufficiently  smaller 
than  and  the  basic  matrix  is  not  too  badly  conditioned,  iterative 

refinement  can  be  expected  to  improve  the  approximations  to  u,  tt, 
and  y  essentially  to  full  normal-precision  accuracy.  In  other  terms, 

T],  01,  and  cp  can  each  be  reduced  to  a  small  integer  multiple  of 

by  the  iterative  refinement.  Since  the  given  data  A,  b,  arid 
c  are  only  assumed  good  to  normal-precision  accuracy,  there  is  no 
reason  to  produce  approximations  to  u,  n,  and  y  having  more  accuracy 
than  that  provided  by  iterative  refinement.  If  difficulties  arise  which 
cannot  be  cleared  up  by  iterative  refinement,  it  is  reasonable  to  argue 
that  they  should  remain  unresolved  pending  acquisition  of  more  accurate 
data. 

We  summarize  these  remarks  by  presenting  a  variation  on  the  simplex- 
method  cycle*. 

Main  Section 

1.  Produce  u  +  6u  and  T]  • 

If  +  6u^  <  T)  for  any  i,  restore  the  basis  used  for  the 
previous  execution  of  the  cycle,  and  go  to  step  8. 


I 

I 

I 

1 

.1 

.! 

J 


2. 

Produce 

TT  + 

6n 

and 

03  . 

3. 

Produce 

z  + 

6z 

and 

5  • 

Let  z  + 

6z 

and 

1 

be  the  corresponding  values  produced 

during  the  preceding  execution  of  the  cycle.  If 
z  +  6z  +  5  <  z  +  6z  -  restore  the  preceding  basis,  and  go 
to  step  8. 

4.  Produce  c.  +  6c.  and  £.  for  each  j  l  {v  .....  v  }, 

J  .1  J  1  m 

1  <  j  <  n  .  For  one  such  value  of  j  set 


t  =  c .  + 
J 


6c 


•8 


and  s  -  j  . 


:i 


For  each  other  such  value  of  j  in  turn  set 


t  =  c .  +  6c.  and  s  -  j 

j  J  J 


if 


:i 

3 

] 

] 

1 


c  .  +  6c .  +  £  .  <  t 
j  j  : 


Having  finished,  if' 


c  +  6c  >  -  £  , 

s  s  —  s 


then  go  to  step  9* 

5.  Produce  y  +  6y  and  cp  . 

If  yi  +  6y^  <  cp  for  all  i,  go  to  step  9- 
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6.  For  one  value  of  i  such  that  +  6y^  >  cp  let 


u  +flu  -T1 

t  =  — Ti — t~  and  r  =  i  . 

yi+oyi'+xp 


For  each  other  such  value  of  i  in  turn  set  t  and  r  as 


above  whenever 


u^u 


y^y^cp 


<  t  . 


7.  Drop  A  from  the  basis  and  add  A  .  Then  form  a  corre- 
v  s 

r 

sponding  basic  matrix,  update  a  decomposition,  or  do  whatever 
else  is  necessary  to  reflect  the  basis  change,  and  go  to 


step  1. 


Refinement  Section 


8-  Form  the  necessary  basic  matrix,  compute  a  decomposition,  and 

produce  u  +  6u,  tt  +  6tt,  y  +  6y  . 

9.  Refine  u  +  6u  and  produce  a  corresponding  error  bound  T|  • 

If  u.  +  fiu,  <  T)  for  any  j,  computation  should  be  suspended, 

J  J 

and  an  indication  should  be  made  that  the  problem  seems  to 
be  degenerate. 

10.  Refine  tt  +  6rr  and  produce  a  corresponding  error  bound  co  . 

11.  For  each  j  j.  •  •  •  >  vm) ,  1  <  i  <  u  produce  c.  +  6c^ 

and  C .  from  the  refined  tt  +  6tt  and  o>  .  Find  s  as  in 

step  4.  If  "c  +  6c  >  C  ,  the  optimal  basis  has  been  found, 
s  s  —  s 
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If  C  >  c  +  6c  >  -  C  ,  the  current  basis  is  prtbabLy, 
though  not  necessarily,  optimal. 

12.  Refine  y  +  6y  and  produce  a  corresponding  error  bound  cp  . 

If  yi  +  6y^  <  tp  for  al.L  i,  computation  should  be  suspended, 
and  an  indication  should  be  made  that  the  objective  function 
appears  to  be  unbounded  from  above.  Otherwise  go  to  step  6. 

We  should  note  in  closing  that,  if  an  implementation  of  the  simplex 
method  constructed  along  the  lines  given  above  produces  a  vector  x  +  fix' 
as  the  computed  approximation  to  an  optimal  solution,  there  is  no  way 
of  checking  its  closeness  to  any  true  optimum  for  the  problem  using 
the  information  at  hand.  This  is  due  to  the  fact  that  the  components 
of  optimal  solutions  to  linear  programming  problems  need  not  vary  con¬ 
tinuously  with  perturbations  in  the  data.  And  round-off  errors  committed 
during  computation  of  the  basic  feasible  solutions  may  be  expressed,  we 
have  seen,  as  perturbations  of  the  basic  matrices,  and  hence  of  the  data. 
A  bound  on  6~  is  useful  only  for  delimiting  the  distance  from  some 
basic  feasible  solution  x"  and  the  computed  vector  x  +  6x"  .  The 
problem  can  quite  possibly  have  a  unique  optimum  x*  for  which  the 
distance 

UxMxwoa 

is  greater  than  the  bound  on  j[6x||  .  Such  a  situation  might  arise,  for 
example,  in  a  problem  having  one  optimal  point  and  one  or  more  widely 
separated  nearly  optimal  points.  Then  round-off  errors  can  quite 


possibly  cause  the  computed  approximation  to  one  ol‘  the  nearly  opt  mu  1 
points  to  have  a  higher  associated  value  ol’  the  objective  function  than 
does  the  approximation  to  the  optimal  point. 


It  is  our  final  goal  to  compare  the  computational  effort  required 
by  each  of  che  three  simplex-method  implementations  which  have  been 
studied  These  implementations  differ  only  in  the  manner  in  wnich  each 
produces  the  vectors  u,  it,  and  y  .  Hence,  a  measure  of  the  worK 
done  by  each  to  produce  these  three  vectors  will  constitute  a  satisfac¬ 
tory  comparison. 

Traditionally  the  number  of  multiplications  and  divisions  (prciuct 
operations ;  needed  to  carry  out  a  computation  has  been  the  measure  of 
the  work  needed  for  that  computation.  The  reason  fcr  this  lies  in  .he 
much  slower  speed  at  which  earlier  generations  of  computers  have  per¬ 
formed  product  operations  in  contrast  to  their  speed  of  performing 
additions  and  subtractions  (.summation  operations).  As  machines  are 
built  with  more  nearly  uniform  speeds  for  all  arithmetic  operations, 
this  tradition  is  becoming  indefensible.  In  the  programs  discussed 
here,  however,  computations  are  composed  of  inner  products  or  of 
additions  of  a  multiple  of  one  vector  to  another-  Thus,  each  product 
operation  generally  gets  paired  with  a  summation  operation,  and  the 
convenience  of  counting  only  the  product  operations  may  be  kept. 

We  begin  by  considering  the  standard  LU  implementation,  since 
several  of  the  operation  counts  which  we  derive  for  it  are  applicable 
tc  the  modified  LU  implementation  as.  well' 

Consider  solving  a  lower  triangular  system 

lv  -  g 
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|$' '  TUI 


if  the  first  k  -  1  components  of  v  are  known.  The  components 

v.  , ...,  v  are  produced  by  setting 

k  m 

vi  *  gi  ‘  'u’j 

in  order  for  i  =  k , . . .,  m  .  This  requires 

(k-i)  +  m  ♦...+  (m-i)  .  "(-iHfc-iHfc-.a? 

product  operations. 

Consider  solving  a  lower  triangular  system 

T 

U  v  =  g 

if  the  first  k  -  1  components  of  v  are  known.  Here,  for  i  =  k,. 

in  order  we  set 


v. 

i 


u.  .V . )/u. . 
1J  j"  11 


) 


so  that  the  cost  to  produce  v  is 


(kn) 

product  operations. 

To  produce  the  vector  u  in  the  standard  LU  implementation,  we 
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solve,  in  order, 


Lt  ~  b  , 


>**»««*  II 


To  produce  the  vector  tt  we  solve,  in  order, 


U  w  =  v 


L  T7  =  W  • 


But  the  matrix  IT  is  identical  in  its  principal  submatrix  of  order 
Tq  -  1  with  the  corresponding  matrix  used  on  the  preceding  execution 
of  the  simplex-method  cycle.  Similarly,  the  vector  y  does  not  differ 
in  its  first  rQ  -  1  components  from  the  corresponding  vector  used  on 
the  preceding  cycle  .  Thus,  as  was  done  in  producing  u,  an  old  copy 
of  w  can  be  used  to  cut  the  cost  of  computing  tt  to 


m(^l)-r0(r0-l)  ^  m(m-1) 


product  operations. 

No  saving  in  the  production  of  y  can  be  arranged.  The  total 
cost  required  for  the  computation  will  be 


m(m-l)  ,  ra(m+lj 
2  2 


product  operations. 

Finally,  if  the  execution  of  the  simplex-method  cycle  results  in 
the  change  of  column  r^  of  the  basic  matrix,  we  must  count  the  cost 
of  updating  the  LU  decomposition.  The  areas  which  change  are  shown 
shaded  next: 


i 


'/Iff 


Now,  a  partial  column  ol‘  the  LU  decomposition, 


is  computed  by  first  letting  i  -  p,,  y,  +  1,...,  v  in  order  and  setting 


i-1 

i.  =  b .  -  T'  t .  .  u . 

iv  iv  M  ij  J\ 


which  requires 


z  ,  .  ,  v  ( v-l)- (u-l)  (li*2 ) 

(|i-lj  +  u  +•••  +  (v-l)  =  — * - -g** — ^ - *- 


I 


product  operations,  and  then  let  tin/*.  *  r  v  '*  L , . .  -  ■  m  n  o-dcr  and 


setting 


l  ~  (b.  -  y  t  . u .  )/u 

iv  v  iv  A  i.)  Jv"  vv 

ij 


which  requires 


(m-v)v 


product  operations. 

To  completely  update  the  LJ  decomposition  oi  the  basic  matrix, 
it  is  required  that  column  r^  be  fully  recomputed,  and  that  ccLimns 
+  1  through  m  be  recomputed  from  row  number  r^  on  down.  This 
costs 


•i(ri  *l) 

-i-gi -  +  (m-r1)r1 


m  k(k-l)-(r,-l)(r  -2) 

£  - - i—  +  (m-k)k' 


k=ri+1 


product  operations. 

Using  the  equalities 


£  k(H~U  ^  |(m+l)m(m-i)-(ri+l)r1(r1-:. )' 


k=r  +1 


I 


m  mm,. 

£  (m-k)k  =  m  £  k  -  £  k‘ 


k-r1+l 


k=r^+l  ksr^+1 


m  m 


=  (m-l)  V  k  -  V  k(k-l) 
k»Vl  kvrxM 


=  (m(m+l)-r1(r1+l)]  -  ^  [m(m+l) (m-l)-r1(r1+l) (r^ -1 )] 


we  find  that  the  cost  of  updating  is 


r,  (r  -.1.) 

- ^ -  +  r](m-r1)  +  £  [  (m+l)m(m-l)-(r1+l)r;L(r1-l)] 


(r  -2)(r1 -l) (m-r  )  , 

_i - £ - —  +  ~  [m(m+l)-r;i(r1+l)l 

^  [(m+l)m(m-l)-(r1+l)r1(r1-l)] 


=  T  '  mri  +  !  ri  +  2mri  "  1  ri  '  5  m  +  I  ri 


product  operations.  Thus,  to  fully  carry  out  the  simplex-method  cyc]e 
using  the  standard  LU  implementation  requires 


2  3  „  5  2  4  5 

-  “l  +  5  rl  +  2mrl  -  2  rl  -  5  m  5  rl 


nr  2 
3 


m(m-l)-(r0-l)(r0-2)  m(m+l) 
+  2  2 
m(mn)-r0(r0-l)  m(m.1) 

2  2 

m(m-l)  +  m(m+l) 

2  2 
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product  operations,  which  simplifies  to 


m  2  .  2  3  .  _  5  2 

T  '  1  4  3  rl  4  2mri  •  2  rl 


4  .5 

i'Vl 


+  5»>2  -  r20  *  2r0 


product  operations. 

This  expression  is  not  as  useful  as  we  would  like  for  comparison 
purposes,  depending  as  it  does  upon  the  present  and  preceding  exchange 
through  the  numbers  r^  and  .  However,  it  is  reasonable  to  assume 
that  each  value  for  r^  and  r^  between  1  and  m  is  equally  likely 
to  occur.  Hence,  the  expected  number  of  product  operations  in  any  execu¬ 
tion  of  the  simplex-method  cycle  would  be 


v2+(2m  +  ^2) 


J 


or 


3 

T  4 


3m 


m 

•  z 


product  operations 

For  the  modified  LU  implementation  some  savings  can  be  made  in 
the  computation  of  the  vectors  u  and  tt  like  the  savings  made  in 
the  standard  LU  implementation. 
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To  produce  the  vector  u  we  solve 


I 


I 


I 


I 


i  - 


3 


3 

3 

1 

3 

I 

] 

1 

I 


Lt  -  b 


and  then  set 


-  irlk]n^k]  ...  r^n^]  ...  Lr'!?n;L|  ... 

m-l  m-I  r.  r,  m~t  m-l  r,  r 


1) 


k  k 


1.  rl 


It 


and  lastly  solve 


,.(k)  (k) 

IJ  u  r.  q 


However,  on  the  previous  execution  of  the  simplex -method  cycle  we  ha 
already  solved 


Lt  ~  b  , 


and  have  set 


( k-l_)  a  [r(k-l)n(k-l)  _  (k-l)n(k-l)j  fr(L)n(l)  (i) 

m-l  m-l  r,  ,  r,  ,  m-l  m-l  r 


U) 


k-1  k-1 


,(k-l) 


Hence,  if  q  is  saved,  we  need  only  set 


q(k)  B  [r(RVkJ  ...  r^kVk^  jq^k_1';  , 

m-l  m-l  r.  r.  ’ 


k  k 


n: 


. — * —  —  . . . 


and  then  solve 


I 


„(k)  00 

U  7u  -  q 


(k) 

The  production  of  q'  costs  m  -  r  product  operations,  and  the 

K 

production  of  u  costs  m(m+l)/2  . 

To  compute  rt ,  we  solve 


w  =  y  , 


I 


then  we  set 


...  ]  ... 

y  1  r,  r,  m-1  m-1  r,  r. 


1  ‘1 


k  k 


•  n lw 


and  solve 


tT  (k) 

L  n  =  p 


Noting,  as  with  the  standard  LU  decomposition,  what  things  remain 
unchanged  from  one  cycle  to  the  next,  we  find  that  the  first  step 
need  cost  only 


m(m+l)-rk(rk-l) 
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I  ! 


t 


I 

I 


i 


i...  ^  . . -  . 


LaWMUlHMli 


in  order,  we  would  carry  out  the  updating  by: 


1.  interchanging  rows  a  and  a  +  1  it'  necessary, 


2.  making  one  division  to  compute  the  Gauss  multiplier 

3.  subtracting  g^*^  times  the  a*'*"1  row  of 

o 

,  ,th 

a  +  1  row. 

In  making  this  last  step,  only  components  a  +  1  through  m 

St 

o+l  row  need  to  be  computed.  Thus,  for  each  value  of  o 
m  +  1  -  o  product  operations  for  a  total  of 


(k+1) 

cr  ' 

o 

from  the 


> 


of  the 
we  require 


(m+l-r^+1)  +...+  (m+l-m+l)  =  2  +  3  +  •••  +  (m+l-rk+^) 

(m-rk+1)(m-rk+1+3) 

2 


Thus,  the  measure  of  the  work  needed  to  carry  out  the  kth  execu¬ 
tion  of  the  simplex-method  cycle  using  the  modified  LU  implementation 
is 


2 

m 


T 


mr 


k+1 


+ 


2 

m 

2 


+  m  - 


r 


k 


+  km 


+  km 


m 

2 
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product  operations  on  an  average  execution  of  the  simplex-method  cycle 
using  the  modified  LU  decomposition. 

For  the  standard,  explicit-inverse  implementation  of  the  simplex- 
method  cycle  described  in  chapter  3,  the  order  of  proceeding  is  slightly 

i 

different.  The  cycle  begins  with  u  and  tt  already  at  hand,  and 
the  vector  y  is  produced  by  forming  the  matrix- vector  product 


B-\  >  y  , 


which  costs  m  product  operations. 

After  the  exchange  of  column  A  for  A  in  the  basis,  columns 

S  V 

r 

1  through  m  of  the  matrix  D  must  be  updated  by  applying  to  each 


she  Gauss- Jordan  elimination  which  reduces  the  vector 


to  the  rth  unit  vector.  That  is,  the  rt^-  row  of  D  is  divided  by 

J" 

tf  ;  then,  for  i  =0,...,  r  -  1,  r  ♦  1,  . m,  t  times  t.he  r  row 

t  h 

is  subtracted  from  the  i  row  of  D  .  In  alL  the  updating  requires 
m(m+-l)  product  operations. 

Finally,  the  same  Gauss-Jordan  elimination  is  applied  to  the 
vector 


at  a  cost  of  m  +  1  product  operations. 

Thus,  one  execution  of  the  s implex -method  cycle  using  this  explicit 
inverse  implementation  costs 


+  2m  +  1 


product  operations.  This  is  to  be  compared  with  the 

IE?+  5jU  km  -  k  -  1 

+■ 

product  operations  needed  on  the  average  to  perform  the  k'’  execution 
of  the  simplex -method  cycle  using  the  modified  LU  implementation.  And 
it  is  to  be  compared  with  the 
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1 


Ell 

C 


4  3  m2 


m 

C 


product  operations  needed  on  tne  average  to  carry  out  the  s imp I  ox -met  hod 
cycle  using  the  standard  LU  implementation. 

Since  k  does  not  generally  grow  larger  than  a  small  multiple  of 
m,  the  modified  LU  implementation  does  not  require  significant  ly  mom 
work  than  the  classical  implementation.  Guaranteed  numerical  stuhi.ity, 
in  short ,  needn’t  be  expensive. 
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